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Abstract 

This  thesis  applies  modern,  multi-variable  control  design 
techniques,  via  a  FORTRAN  computer  algorithm,  to  U.S.  Army 
helicopter  models  in  hovering  flight  conditions. 
Eigenstructure  assignment  and  Linear  Quadratic  Regulator 
(LQR)  theory  are  utilized  in  an  attempt  to  achieve  enhanced 
closed  loop  performance  and  stability  characteristics  with 
full  state  feedback.  The  addition  of  cross  coupling  weights 
to  the  standard  LQR  performance  index  is  specifically 
addressed.  A  desired  eigenstructure  is  chosen  with  a  goal 
of  reduced  pilot  workload  via  performance  characteristics 
and  modal  decoupling  consistent  with  current  helicopter 
handling  qualities  requirements.  Cross  coupling  weighting 
is  shown  to  provide  greater  flexibility  in  achieving  a 
desired  closed  loop  eigenstructure.  Also,  while  the 
addition  of  cross  coupling  weighting  is  shown  to  eliminate 
stability  margin  guarantees  associated  with  LQR  methods,  the 
study  shows  that  the  modified  algorithm  can  achieve  a  closer 
match  to  a  desired  eigenstructure  than  previous  versions  of 
the  program  while  maintaining  acceptable  stability 
character i st i cs . 
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HELICOPTER  FLIGHT  CONTROL  SYSTEM  DESIGN  USING  THE  LINEAR 
QUADRATIC  REGULATOR  FOR  ROBUST  EIGENSTRUCTURE  ASSIGNMENT 

I.  Introduction 

This  thesis  applies  modern,  multi-variable  control 
design  techniques,  via  a  FORTRAN  computer  algorithm,  to  U.S. 
Army  helicopter  models  in  hovering  flight  conditions. 
Eigenstructure  assignment  and  Linear  Quadratic  Regulator 
(LQR)  theory  are  utilized  in  an  attempt  to  achieve  enhanced 
closed  loop  performance  and  stability  characteristics  with 
full  state  feedback.  The  addition  of  cross  coupling  weights 
to  the  standard  LQR  performance  index  is  specifically 
addressed.  A  desired  eigenstructure  is  chosen  with  a  goal 
of  reduced  pilot  workload  via  performance  characteristics 
and  modal  decoupling  consistent  with  current  helicopter 
handling  qualities  requirements.  This  work  builds  upon 
previous  efforts  by  Robinson  [1]  and  Huckabone  [2]  at  the 
Air  Force  Institute  of  Technology  (AFIT),  the  major 
difference  being  the  addition  of  cross-coupling  weighting  in 
the  LQR  performance  index. 

He1.lcop.tgrs 

Helicopters  serve  as  a  good  platform  for  the 
application  of  automatic  control.  The  aerodynamic  and 
structural  design  of  rotary  winged  aircraft  create  inherent 
problems  with  regard  to  stability  and  control.  In  order  to 
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apply  the  tools  of  control  design,  a  mathematical  model  of 
the  system  is  needed.  McRuer  [3: chap  4]  presents  the 
assumptions  and  techniques  used  to  represent  an  aircraft  as 
a  set  of  linear  differential  equations  with  constant 
coefficients,  or  equations  of  motion.  Reid  [4:chap  6] 
provides  the  techniques  by  which  these  equations  can  be 
represented  as  transfer  functions  or  via  the  state  space 
model 


X  =  Ax  *  Bu 


(1) 


The  state  space  model  is  used  for  the  design  technique 
presented  in  this  thesis. 

Two  helicopter  models  are  used  in  this  thesis.  One  is 
based  on  the  AH-64  (Apache)  attack  helicopter  [5]  the  other 
is  the  UH-60  (Blackhawk)  utility  helicopter  [6].  Both 
aircraft  are  considered  to  be  conventional  helicopters  in 
that  they  have  a  single  main  rotor  to  produce  primary  lift 
and  a  single  tail  rotor  to  produce  anti -torque  and 
directional  control  forces.  The  models  used  here  represent 
the  aircraft  in  hovering  flight. 

The  state  vector  x  in  equation  (1)  is  the  column  vector 
comprised  of  the  following  respective  states: 

u  -  longitudinal  (long)  velocity  (vel)  [feet/second] 

V  -  lateral  (lat)  velocity  [feet/second] 
w  -  vertical  velocity  [feet/second] 
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p  -  roll  rate  [radians/second] 
q  -  pitch  rate  [radians/second] 
r  -  yaw  rate  [radians/second] 

<t>  -  roll  angle  [radians] 

6  -  pitch  angle  [radians] 

The  input  vector  u  is  the  column  vector  comprised  of  the 
following  respective  inputs: 

6e  -  collective 
6b  -  longitudinal  cyclic 
6,  -  lateral  cyclic 
6^  -  tail  rotor 

Collective  input  refers  to  a  simultaneous  change  in  the 
angle  of  attack  of  all  the  main  rotor  blades  that  causes  a 
change  in  magnitude  of  the  lift  produced  by  the  main  rotor. 
Cyclic  input  refers  to  independent  angle  of  attack  changes 
of  the  blades  that  result  in  a  change  of  main  rotor  lift 
direction.  Tail  rotor  input  is  similar  to  collective  input 
applied  to  the  tail  rotor.  For  the  AH-64  model,  the  inputs 
are  expressed  in  degrees  of  blade  movement,  while  the  UH-60 
inputs  are  expressed  as  inches  of  control  movement. 

Standard  sign  conventions  for  states  and  controls  are  used 
as  presented  in  reference  [3]. 

The  A  and  B  matrices  of  equation  (1)  contain  the 
constants  of  the  equations  of  motion.  These  matrices  are 
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presented  in  appendices  A  and  B  for  the  AH-64  and  UH-60, 
respectively.  The  AH-64  model  has  been  non-dimensional i zed 
as  described  by  Huckabone  [2:52-54]. 

The  purpose  of  designing  an  automatic  control  system 
for  an  aircraft  is  to  enhance  its  handling  or  flying 
qual ities. 

Flying  qualities  determine  the  ease  and  accuracy 
with  which  a  pilot  can  accomplish  the  various 
tasks  or  maneuvers  that  constitute  the  aircraft's 
mission.  The  elements  which  directly  affect 
flying  qualities  are  the  stability  and  control 
characteristics  of  the  aircraft  which  link  the 
controllers  that  the  pilot  manipulates  to  the 
aircraft  response  states  that  the  pilot  desires  to 
control  [7:2-1]. 

Factors  that  can  cause  adverse  flying  qualities  are: 
instabilities,  sluggish  response  and  uncommanded  responses. 
Specific  flying  quality  requirements  for  helicopters  are  set 
forth  in  Aeronautical  Design  Standard-33  (ADS-33)  [8]. 
Analysis  of  the  open  loop  eigenstructure  of  the  helicopter 
model  will  reveal  the  presence  of  instabilities  and  state  to 
state  coupling.  The  B  matrix  of  the  helicopter  model  shows 
the  control -state  coupling  that  is  specifically  addressed  in 
this  thesis. 

The  need  to  improve  handling  qualities  arises  from  the 
desire  to  reduce  the  pilot's  workload  as  much  as  possible. 
This  is  done  to  allow  the  pilot  to  perform  tasks,  other  than 
flying,  in  conjunction  with  the  aircraft's  mission. 

Examples  of  tasks  include  navigation,  communication  and 
weapon  delivery.  One  way  that  flying  qualities  may  be 
improved  is  by  the  reduction  of  control -state  coupling  or 
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cross  coupling.  Ideally,  a  helicopter  pilot  should  have  the 
ability  to  command  precise  control  of  the  hovering 
aircraft's  position  via  direct,  single  axis  control.  The 
desired  control -state  coupling  Is: 


Collective 

Longitudinal  Cyclic  - 
Lateral  Cyclic 
Tall  Rotor 


Vertical  Velocity 
Longitudinal  Velocity 
Lateral  Velocity 
Yaw  Rate 


Automatic  Flight  Control  Systems  (AFCS)  are  currently 
utilized  In  helicopters  to  reduce  pilot  workload.  One  of 
the  most  modern  aircraft  designs,  the  MH-60K,  employs  the 
following  features  to  augment  flying  qualities  In  a  hover: 

1.  Pitch,  roll  and  yaw  stability 

2.  Cyclic,  collective  and  directional  trim 

3.  Pitch  and  roll  attitude  hold 

4.  Heading  hold 

5.  Altitude  hold 

6.  Coupled  hover  (Inputs  through  automatic 
pilot)  [9] 

These  functions  greatly  enhance  the  mission  effectiveness  of 
the  aircraft  and  crew  In  a  combat  aircraft  employing  many 
complicated  weapons  systems. 
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Flight  Control  Systems  Design 


AFCS  for  helicopters  are  primarily  designed  using 
classical  or  single  input,  single  output  (SISO)  methods. 
These  methods  require  that  the  multiple  input,  multiple 
output  (MIMO)  state  space  model  of  the  system  be  broken  down 
into  scalar  subsystems  or  transfer  functions.  Even  recent 
research  by  Osder  and  Caldwell  concludes  that  “A  practical 
process  for  designing  such  multiple  input,  multiple  output 
helicopter  systems  starts  by  decoupling  controls  into  four 
single  input,  single  output  axes  ...  [10].“  The 
simplification  of  the  MIMO  system  is  deemed  necessary 
because  SISO  techniques  become  very  complicated  when  all 
input-output  relationships  are  addressed,  especially  cross 
coupling  of  controls  and  states.  This  normally  leads 
control  system  designers  to  limit  their  focus  to  the  worst 
cross  coupling  areas  while  ignoring  the  rest. 

As  an  example,  the  UH-60  helicopter  flight  control 
system  incorporates  a  mechanical  mixing  unit  designed  to 
eliminate  coupled  control  inputs  [11].  The  mixing  unit 
links  only  four  out  of  the  twelve  possible  cross  couplings 
of  the  desired  control -state  matches.  This  is  done  even 
though  coupling  is  present  in  all  of  the  control -state 
pairings. 

Modern  design  techniques  allow  for  direct  use  of  the 
MIMO  state  space  model  in  control  system  design  algorithms. 
There  is  no  need  to  break  down  an  integral  model  into  SISO 
subsystems.  Thus  all  the  control  state  relationships. 
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including  cross  couplings,  are  considered  in  the  design 
process.  Two  popular  MIMO  design  techniques  are 
eigenstructure  assignment  and  the  Linear  Quadratic 
Regulator. 

Eigenstructure  assignment  allows  control  system 
designers  to  prescribe  desired  closed  loop  eigenvalues  and 
eigenvectors  for  a  given  system,  thus  achieving  desired 
performance  characteristics.  Garrard  and  Liebst  used 
eigenstructure  assignment  "to  design  a  feedback  system  for 
use  in  precise  hovering  control  for  a  modern  attack 
helicopter.  Eigenvalue  placement  is  used  for  stability 
enhancement  and  eigenvector  shaping  is  used  for  modal 
decoupling  [12]."  But,  as  Moore  has  shown,  eigenstructure 
assignment  does  not  provide  a  unique  control  system  design 
[13].  This  flexibility  allows  the  designer  to  augment 
eigenstructure  assignment  with  additional  design  methods. 

The  Linear  Quadratic  Regulator  (LQR)  is  an  optimal 
control  design  method  that  yields  a  full  state  feedback 
controller  which  minimizes  the  quadratic  performance  index 

m 

♦T-J  (x^Oic  +  u^Mu) dt  (2) 

0 


where  Q  is  the  symmetric,  positive  semi -definite  state 
weighting  matrix  and  R  is  the  symmetric,  positive  definite 
control  weighting  matrix.  Note  that  this  performance  index 
is  quadratic  in  both  state  and  control  variations  from 
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nominal  conditions,  hence  minimizing  attempts  to  maintain 
small  plant  deviations  with  small  control  inputs.  The  LQR 
method  can  provide  a  robust  closed  loop  solution  with 
guaranteed  stability  margins  [14:sect  6]. 

Using  the  eigenstructure  assignment  method,  the  control 
system  designer  can  specify  the  desired  performance  of  a 
given  system.  The  LQR  method  provides  a  robust  design 
solution.  Combined  use  of  these  methods  yields  the  best  of 
both  worlds. 

Program  Background 

Captain  Jeffrey  0.  Robinson  developed  an  algorithm  at 
AFIT  which  utilized  eigenstructure  assignment  and  the  LQR 
method  in  defining  linear,  full  state  feedback  gain  for 
aerospace  systems  [1].  Robinson's  algorithm  was  written 
exclusively  for  use  with  the  MathWorks  Inc.  software  package 
MATLAB^"  [15]  and  was  limited  to  eigenvalue  assignment 
only.  Captain  Thomas  C.  Huckabone,  also  at  AFIT,  augmented 
Robinson's  work  by  introducing  eigenvector  assignment  [2]. 
Huckabone  rewrote  Robinson's  program  in  FORTRAN  using 
MATLAB^"  only  as  a  method  of  manipulating  input  and  output. 

In  addition  to  newly  written  routines,  Huckabone 's  work 
utilized  existing  subroutines  from  the  LQGLIB  [16]  and  IMSL 
[17]  packages  available  on  the  AFIT  computer  system. 
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The  algorithm  minimizes  the  performance  Index 


where 

n 

Fe 

X, 

X. 

Vd 

V. 

Fv 

e 


number  of  states 
eigenvalue  weighting  matrix 
desired  closed  loop  eigenvalue 
achieved  dosed  loop  eigenvalue 
desired  closed  loop  eigenvector 
achieved  closed  loop  eigenvector 
eigenvector  weighting  matrix 
eigenvector  minimization  constant 


The  feedback  controller  Is  obtained  via  LQR  methods,  which 
minimize  the  LQR  performance  Index  presented  In  equation 
(2).  Specifically,  the  algorithm  varies  the  LQR  performance 
Index  state  and  control  weighting  matrices,  Q  and  R 
respectively,  using  an  optimization  method  based  upon  the 
Nelder-Mead  simplex  algorithm  presented  In  reference  [18]. 

This  thesis  Is  a  direct  extension  of  Huckabone's  work. 
His  algorithm  is  modified  so  that  the  cross  coupling 
weighting  matrix  S  Is  Included  In  the  LQR  performance  Index. 


9 


The  performance  index  is  now  written  as 


(4) 


The  modification  was  added  to  provide  increased  flexibility 
in  achieving  the  desired  eigenstructure.  The  algorithm  is 
then  applied  to  mathematical  models  of  conventional,  modern 
helicopters  in  hovering  flight  conditions.  The 
eigenstructure  and  stability  characteristics  achieved  via 
the  modified  algorithm  are  compared  with  results  that  do  not 
utilize  the  cross  coupling  weighting  matrix. 
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u.  Cross  CouDlInq  Extension 


Theory 

The  majority  of  the  mathematical  theory  and  equations 
necessary  for  development  of  the  algorithm  are  reported  In 
detail  by  Huckabone  [2:9-22].  In  particular,  Garrard, 

Llebst  [12]  and  Moore  [13]  provide  elgenstructure  assignment 
theory  while  Ridge 1y  [14]  provides  LQR  theory.  Some  general 
equations  are  presented  below  for  convenience.  The  theory 
necessary  to  Introduce  the  cross-coupling  weighting  matrix  S 
Is  presented  here  In  detail. 

General .  Again,  the  standard  state  equation  of  a 
multivariable,  linear,  time  Invariant,  feedback  system  Is 

X  =  Ax  *  Mu  (5) 


Assuming  full  state  feedback  and  a  B  matrix  with  full  column 
rank  yields  a  linear,  feedback  control  law  of 

n  =  -JOT  (6) 

Again,  the  standard  LQR  method  Involves  minimizing  the  cost 
function 


«7- J  (jr'Or  ♦  a  'Jte)  dt 

0 


(7) 
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where  Q  is  the  symmetric,  positive  semi-definite  state 
weighting  matrix  and  R  is  the  symmetric,  positive  definite 
control  weighting  matrix.  From  LQR  theory,  an  optimal 
feedback  gain  matrix 

(8) 


is  obtained  where  the  symmetric  matrix  P  is  the  stabilizing 
solution  to  the  algebraic  Ricatti  equation 

PA  +  A*P  +  0  -  Pn-^B*P  -  0  (9) 


Cross-CouDi ino.  Introduction  of  the  performance  index 


(10) 


allows  for  weighting  and,  consequently,  minimization  of 
combined  state  and  control  terms  (i.e.  x,Uj)  via  the 
matrix  S.  The  original  LQR  performance  index  in  equation  (7) 
allows  minimization  of  pure  state  and  control  terms  (x,Xj  or 
u,Uj)  only^.  The  integrand  of  equation  (10)  can  be  expanded 
via  matrix  multiplication  to 

x'gr  +  *  u*Pu  (11) 


^The  standard  LQR  performance  index  is  equivalent  to  the  new 
performance  index  when  all  elements  of  the  matrix  S  are  zero, 
hereafter  referred  to  as  S=[0]. 
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Slnpllfylng  equation  (11)  via  scalar  addition  and 
substituting  back  Into  equation  (10)  yields  the  rewritten 
performance  1 ndex 

J  (x'l^  +  2jr1ita  n  *Mu)  dt  ( 12 ) 

0 


Anderson  and  Moore  [19:56]  show  that  standard  LQR  methods 
can  be  utilized  to  find  an  optimal  feedback  solution  that 
minimizes  a  performance  index  written  in  the  form  of 
equation  (12)  with  the  following  constraints* 


Jl  >  0 

(13) 

(14) 

This  is  shown  by  rewriting  the  Integrand  of  equation  (12)  as 

x*{0-dM^S*)x  *  (15) 


where 


|i  ■  a 


(16) 


Via  substitution,  the  original  state  equation  (5)  is 


*The  expressions  >  0  and  Z  0,  when  used  with  matrices, 
refer  to  the  matrix  being  positive  definite  and  positive  semi- 
definite,  respectively. 
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rewritten  as 


i  =  U  -  Btt-^a*)x  +  J|i 


(17) 


The  redefined  LQR  problem  yields  an  optimal  control  law  of 

|1  =  (18) 

where  P  is  now  the  stabilizing  solution  to  the  algebraic 
Ricatti  equation 

( O-SR-^a*)  =0  ( 19 ) 

The  optimal  control  law  for  the  original  system  is  then 
found  by  combining  equations  (16)  and  (18)  which  results  in 

u  =  (B'p  *  a*)  X  ( 20 ) 


which  satisfies  the  following  necessary  and  sufficient 
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condition  for  a  positive  definite  matrix  [20:331] 

for  all  vectors  X  ^  0  (22) 


Defining 


(23) 


and  substituting  the  partitioned  matrix  of  equation  (21) 
into  the  inequality  (22)  yields 

x^Ok^  +  2s^Sx^  *  XtXx^  >  0  (24) 

Setting  Xt=0  yields 

x,*Mx,  >  0  (25) 

when  Xt^O  ,  which  meets  the  necessary  and  sufficient 
condition  for  positive  definiteness  of  matrix  R,  thus 
ensuring  the  constraint  of  equation  (13)  is  met.  Rewriting 
the  inequality  (24)  as 

^(Q-n^8^x^  +  >  0 


where 

X  •  x^  *  M~*0*x^  (27) 
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and  setting  x=0»  Xi#0  yields 

>0  (28) 

which  meets  the  necessary  and  sufficient  condition  for 
positive  definiteness  of  matrix  [Q-SR~^S^],  thus  satisfying 
the  constraint  of  equation  (14).  A  more  general  proof  for 
satisfying  these  constraints  is  provided  in  Kreindler  and 
Jameson  [21:147]. 

Algorithm  Changes 

The  original  FORTRAN  program  EIGSPACE,  as  written  by 
Huckabone,  was  changed  where  necessary  to  allow  the  use  of 
the  cross  coupling  weighting  matrix  S  as  described  in  the 
theory  presented  above.  In  addition,  the  program  was 
cleaned  up  to  eliminate  unnecessary  procedures.  Details 
presented  in  this  thesis  represent  only  the  major  changes 
made  for  this  thesis.  The  original  EIGSPACE  program  and 
subroutine  descriptions  are  presented  in  detail  as 
appendices  A  and  B  of  Huckabone 's  work  [2:72-103].  The 
LQGLIB  package  of  subroutines  and  the  IMSL  subroutine  were 
not  altered  and  are  described  in  references  [16]  and  [17], 
respectively.  The  modified  program  is  presented  in 
appendix  C.  Operating  instructions  and  options  are 
presented  in  appendix  0. 

The  major  changes  to  the  original  algorithm  occur  in 
two  subroutines.  Subroutine  PP  is  modified  to  allow  use  of 
a  standard  LQR  solver,  in  this  case  the  LQGLIB 
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■  '  '  .  .  .  '■  '  UI-  '  .  .  giJ  Jijiii  Iiii'ii 

subroutine  REG.  Subroutine  MAKEQRS  forms  the  weighting 
matrices  Q,  R  and  S.  Also,  the  SORT  subroutine  was 
eliminated  as  It  was  unnecessary.  Otherwise,  the 
optimization  methods  and  program  flow  for  the  algorithm  are 
unchanged  from  that  presented  In  detail  In  Huckabone's 
thesis  [2:23-30]. 

Using  the  Standard  Regulator  (Subroutine  PP).  As 
stated  In  the  cross  coupling  theory  above,  to  modify  a 
standard  LQR  solver  to  Include  the  cross  coupling  weighting 
matrix  S,  the  Q  and  A  Input  matrices  are  redefined  as 

(29) 


and 


umr  » 


(30) 


respectively.  The  regulator  gain  matrix  of  the  redefined 
problem  Is  then  given  as 

a  -  (31) 

where  P  Is  the  stabilizing  solution  to  the  Ricatti  equation 
(19).  The  optima]  control  law  for  the  desired  problem  is 
given  In  equation  (20)  and  yields  a  regulator  gain  matrix 
for  the  desired  system  of 

JT + S')  (32) 


17 


Thus  the  desired  closed  loop  A  matrix  for  the  original 
system  Is 

MCLB  >  U-BK]  (33) 

Creating  R  and  S  (StibrQMtine  MAKEORS).  The 
parameters  varied  by  the  optimization  subroutines  are 
contained  In  the  vector  XGUESS.  XGUESS  makes  up  the  upper 
triangular  elements  of  the  matrices  QH  and  RN  In  the  case 
where  S  Is  not  varied.  In  the  case  where  S  Is  varied, 

XGUESS  also  contains  all  elements  of  the  matrix  SN.  In  both 
cases,  the  Initial  XGUESS  Is  set  such  that  Q  and  R  are 
appropriately  dimensioned  Identity  matrices  and,  when 
appropriate,  S=:[0]. 

.  In  the  nonvarlable  S  case,  Q  and  R  are  formed  as  In  the 
MAKEQR  subroutine  from  the  original  program  [2:100]  with 
S=[0]  Introduced  as  Input  to  the  PP  subroutine.  The 
following  description  of  how  the  algorithm  creates  the  Q,  R 
and  S  matrices  applies  only  to  the  non-zero  S  case.  While 
setting  SN  equal  to  an  all  zero  matrix  would  properly  form 
the  Q  and  R  matrices  In  the  nonvarlable  S  case,  the 
algorithm  requires  that  the  two  cases  be  handled 
differently.  Specifically,  In  the  non-variable  S  case,  the 
XGUESS  vector  does  not  contain  the  additional  parameters  for 
SN  as  they  would  be  unnecessarily  Iterated  during  the 
optimization  process. 
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As  stated  In  the  natrlx  definiteness  theory,  In  order 
to  allow  a  variable  S  In  the  algorithm,  the  positive 
definite  matrix 


(34) 


must  be  created.  This  matrix,  renamed  QRS  In  the  program, 
can  be  formed  as 

OB0  -  MBV*JMOr  (35) 

where  MHN  Is  the  real  symmetric  matrix 

— ea 

Note  that  equation  (35)  is  equivalent  to 


where  *  denotes  the  Hermit Ian  transpose  of  a  matrix. 

Ridgely  and  Banda  [14:2-7]  provides  that  as  long  as  MHN  Is 
nonsingular,  the  eigenvalues  of  QRS  are  all  positive,  hence 
QRS  is  positive  definite  [20:331].  The  initial  setting  of 
MHN  and  the  subsequent  iteration  process,  virtually  assure 
that  MHN  will  always  be  nonsingular;  however,  the  algorithm 
Is  designed  to  yield  an  error  message  if  this  occurs. 
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Substituting  the  partitioned  NHN  matrix  of  equation 


(36)  into  equation  (35)  yields: 


(38) 


Thus,  the  weighting  matrices  generated  by  the  algorithm  are: 

Q (m*(m  *  aa*8W  (39) 

M  -  *  8a**ai  (40) 

8  -  08*88  *  88*88  (41) 

Eigenvalue  Pairing  (Subroutine  SORT).  As  stated 
previously,  the  algorithm  minimizes  the  performance  index 


where 

n 

Fe 

X. 

X. 

Vd 

V. 

Fv 

e 


number  of  states 
eigenvalue  weighting  matrix 
desired  closed  loop  eigenvalue 
achieved  closed  loop  eigenvalue 
desired  closed  loop  eigenvector 
achieved  closed  loop  eigenvector 
eigenvector  weighting  matrix 
eigenvector  minimization  constant 
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The  above  performance  Index  may  not  accurately  reflect  the 
designer's  desired  value  during  execution  of  the  program. 
This  Is  due  to  the  fact  that  the  algorithm  must  select  which 
desired  and  achieved  eigenvalues  are  to  be  paired  together. 
(Note  that  the  eigenvectors  are  paired  In  accordance  with 
their  respective  eigenvalues.)  Proper  pairing  can  not  be 
guaranteed  by  the  program  as  It  has  no  provisions  for 
Identifying  eigenvalues  by  mode.  Currently,  the  algorithm 
pairs  the  eigenvalues  by  comparing  each  achieved  value  with 
the  value  that  Is  designated  as  the  first  desired  eigenvalue 
and  selecting  the  closest  as  Its  complement.  This  procedure 
continues  for  each  desired  eigenvalue,  according  to  a 
selected  order.  The  original  EIGSPACE  program  attempted  to 
utilize  the  SORT  subroutine  as  a  method  of  selecting  the 
order.  The  subroutine  sorted  the  eigenvalues  In  order  of 
Increasing  magnitude  and  was  used  to  sort  both  the  desired 
and  achieved  eigenvalues. 

Sorting  the  achieved  eigenvalues  was  determined  to  be 
unnecessary  since  the  algorithm  Ignores  the  sorted  order 
when  determining  the  closest  match.  Sorting  the  desired 
eigenvalues  was  also  deemed  unnecessary  as  the  user  may 
select  the  order  when  the  desired  elgenstructure  Is  Input  to 
the  program.  It  Is  Important  to  note  that  the  Input  order 
Is  Important  since  different  orders  can  yield  different 
solutions.  A  simplified  example  Illustrates  this  problem. 

Figure  1  shows  a  possible  mapping  of  achieved  and 
desired  eigenvalues.  It  appears  obvious  that  In  determining 
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how  close  the  achieved  values  are  to  the  desired  structure, 
a  designer  would  pair  the  complex  achieved  eigenvalues  with 
the  complex  desired  eigenvalues  and  the  real  achieved  value 
with  the  real  desired  value.  However,  if  the  real  desired 
value  is  considered  first  in  the  pairing  method  used  by  the 
algorithm,  it  will  pair  off  with  one  of  the  complex  values, 
as  they  are  closer.  Thus  the  performance  index  will  not 
accurately  reflect  what  the  designer  would  like  to  define  as 
closeness.  This  will  affect  the  minimizing  path  that  the 
algorithm  takes  in  searching  for  the  best  solution. 


Figure  1 

Example  of  Eigenvalue  Pairing 
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It  should  also  be  noted  that  the  input  order  of  the 
sign  of  the  imaginary  component  of  desired  complex 
eigenvalues  is  important  as  different  pairings  may  occur  for 
different  orders.  The  algorithm  is  written  so  that  the 
achieved  eigenstructure  determined  by  the  DLQGLIB  subroutine 
EIGVV  always  yields  any  complex  pairs  with  the  positive 
imaginary  eigenvalue  first.  Knowing  this,  the  user  can 
input  the  desired  eigenstructure  accordingly. 

For  the  simplified  example  above,  the  input  order  of 
the  desired  structure  to  avoid  the  mismatching  problem  is 
obvious.  Unfortunately,  higher  state  problems  do  not  always 
provide  the  same  easy  insight  and  the  number  of  permutations 
increases  dramatically.  Thus,  the  input  order  of  the 
desired  eigenstructure  becomes  another  designer  chosen 
parameter. 


23 


ill.  stability  Robustness 


Thflfir.Y 

Ridgely  and  Banda  [14]  present  theory  for  stability 

robustness  of  MIMO  systems  as  well  as  guaranteed  stability 

margins  using  the  LQR  solution.  Specifically,  the  notion  of 

independent  gain  and  phase  margins  is  introduced  as  follows: 

Independent  gain  margins  (IGM)  are  limits  within 
which  the  gains  of  all  feedback  loops  may  vary 
independently  at  the  same  time  without 
destabilizing  the  system,  while  the  phase  angles 
remain  at  their  nominal  values.  Independent  phase 
margins  (IPM)  are  limits  within  which  the  phase 
angles  of  all  loops  may  vary  independently  at  the 
same  time  without  destabilizing  the  system,  while 
gains  remain  at  their  nominal  values  [14:3-73]. 

The  following  relationships  are  shown  to  exist 


^  <  IGM  <  (43) 

+«  l-a 


and 


(I)  <  ^  (t) 


where  a  is  the  minimum  singular  value,  for  all  frequencies, 
of  the  return  difference  matrix  given  by 

a  -  [X+r(i«X-JI) -1*]  (45) 


Note  that  equation  (45)  must  satisfy  the  constraint  a  ^  1. 
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The  Inequalities  (43)  and  (44)  are  conservative,  thus  they 
yield  guaranteed  mlninum  stability  margins. 


Ridgely  and  Banda  [14:sect  7]  also  derive  the  Kalman 
Inequality  from  LQR  relationships,  which  Is: 


(46) 


Guaranteed  LQR  stability  margins  are  derived  from  equation 
(46)  under  the  restriction  R  -  pi,  where  p  Is  any  positive 
scalar.  Hence,  the  Kalman  inequality  (46)  simplifies  to 

iI*K{jaX-A)  i  J  ( 47 ) 

This  can  be  true  only  when 

«  =  ii[X+JC(i«X-A)"‘3l  k  1 


Substituting  a  1  Into  equations  (43)  and  (44)  yields  the 
following  guaranteed  minimum  LQR  stability  margins  under  the 
restriction  R  =  pi. 

1  <  IGH  <  ••  (49) 

2 

-60*  <  IPIf  <  60*  (50) 


'JaH  ii  K-Vi  - 1  □  im*i  ii*!- i  •!  ii.«i 


As  reported  by  Huckabone  [2],  Safonov  and  Athans  [22] 
show  that  when  R  is  diagonal,  the  stability  margins  of 
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equations  (49)  and  (50)  hold  as  long  as  the  perturbations  to 

each  channel  of  the  system  occur  Independently.  Independent 

perturbations  are  Implied  by  an  R  matrix  having  eleieents  of 

the  same  relative  magnitudes.  Huckabone  points  out  that: 

For  the  case  of  any  general  R,  the  Independent 
stability  margins  ...  cannot  be  guaranteed  and 
often  will  go  outside  of  these  bounds.  However, 
as  previously  mentioned,  the  equations  for  IGM  and 
IPM  provide  conservative  values.  While  the  choice 
of  any  general  R  may  not  provide  the  guaranteed 
stability  margins  ...  the  system  may  still  provide 
acceptable  stability  characteristics  [2:22]. 

Introduction  of  the  S  matrix  In  the  algorithm  virtually 

assures  that  R  will  not  be  diagonal.  This  Is  seen  by 

reviewing  the  equation 


ji  >  ^  8a**ai 


(51) 


In  the  algorithm,  the  positive  definite  matrix  RM^RM  can  be 
restricted  via  Input  codes  as  follows: 

r-code  =  1  ^  RM»RM  =  pi 

r-code  =  Z  •*  RM*RH  Is  diagonal 

r-code  =  3  ^  RM*RH  Is  general 

In  the  cases  where  r-code  =  1  or  2,  the  off-diagonal 
elements  of  R  can  only  be  zero  If  SN^*SN  Is  diagonal.  This 
would  require  relationships  between  elements  of  SN.  Since 
each  element  of  SN  Is  Independently  generated,  no  forced 
correlations  between  elements  exist.  While  possible.  It  Is 
highly  unlikely  that  the  elements  of  SN  will  randomly  meet 
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thtt  rsquIrttMints  for  SN^*SII  to  bo  diagonal.  In  tha  case 
whare  r-coda  ^  3,  corral at Ions  batwaan  SN  and  RN  nust  ba 
Mat;  tharefora.  It  Is  evon  more  unlikely  that  R  can  ba 
diagonal . 

Fortunately,  It  Is  possible  that  the  algorithm  will 
provide  an  R  matrix  close  to  diagonal.  Huckabone's  results 
for  the  AH-64  helicopter  yield  an  R  matrix  close  to 
diagonal . 

The  significance  of  this  R  matrix  is  that  because 
it  comes  close  to  approximating  a  diagonal  matrix, 
the  minimum  singular  value  of  the  return 
difference  matrix  is  nearly  one.  It  turns  out 
that  R  being  near  diagonal  is  a  general  result  for 
all  of  the  cases  run  with  this  example  for  R  >  0. 
Therefore  ...  the  resulting  closed  loop  systems 
will  still  possess  good  independent  stability 
margins  [2:59]. 

Again  looking  at  equation  (51),  If  RH*RM  is  close  to 
diagonal,  the  resulting  R  will  be  close  to  diagonal  if 


\  (52) 


While  the  preceding  shows  that  good  stability  margins 
are  possible,  it  also  demonstrates  that  stability  margins 
are  no  longer  guaranteed.  In  fact,  it  can  be  shown  that  the 
eigenvalues  of  a  system  can  be  placed  anywhere  within  the 
stable  region  of  the  complex  plane  using  cross  coupling 
weights. 

Robinson  [1]  demonstrated  how  his  version  of  the 
algorithm  would  find  solutions  only  within  an  achievable  LQR 
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region  for  a  two  state,  single  input  system  defined  by  the 
following  matrices: 

0 

1  , 

Recall  that  Robinson's  algorithm  did  not  use  cross  coupling 
weights  and  matched  eigenvalues  only.  Note  that  for  a 
single  input  system,  R  =  pi  is  satisfied,  thus  the  LQR 
guaranteed  stability  margins  of  equations  (49)  and  (50) 
apply.  In  fact  these  margins  restrict  the  solution,  thus 
preventing  the  algorithm  from  achieving  the  desired 
eigenvalues.  The  algorithm  developed  in  this  thesis,  using 
cross  coupling  weights,  was  able  to  achieve  eigenvalues 
outside  of  that  restricted  region.  In  fact,  using  the 
system  described  by  the  above  A  and  B  matrices,  it  can  be 
demonstrated  that  the  LQR  with  cross  coupling  weights  can 
achieve  any  dosed  loop  stable  solution,  thus  inferring  that 
there  are  no  guaranteed  stability  margins  when  cross 
coupling  weights  are  applied  in  the  algorithm. 

The  closed  loop  characteristic  equation  for  a  system  is 
defined  as: 

det[XX-A4-jyc]  -  0  (53) 
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For  the  system  defined  above  by  A  and  B,  with  K  defined  as: 


JC.[Jc,  Jc,] 


equation  (53)  becomes: 

)L*+J^X+l:j^  -  0 


(54) 


Thus  the  ability  to  arbitrarily  assign  ki  >  0  and  k,  >  0 
would  provide  closed  loop  eigenvalue  assignability  within 
the  entire  stable  (or  left)  half  of  the  complex  plane. 

For  this  example t  let  the  LQR  weighting  matrices  be 
defined  as: 


0 


<Jii  0 
0 


Jl  =  1 


0  . 


The  regulator  gain  and  Ricatti  equations  using  cross 
coupling  weights  are  repeated  here  for  convenience: 

K  ^  +  8*)  (55) 

■)  =0  (  56 ) 

where  P  for  this  example  can  be  defined  as  follows: 

^  [  Pu  Pia 

Solving  for  the  regulator  gain  In  equation  (55)  for  this 
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example  yields: 


*1  -  Pi2  ♦ 

(57) 

^  “  P22 

(58) 

Expanding  the  Ricatti  equation  (56)  yields  the 

relationships: 

following 

^Pl2  *  »l)*  -  Qxi 

(59) 

-  <^22  *  2Pi2 

(60) 

Combining  these  results  yields  the  following  equations  for 

the  scalar  gains  in  terms  of  the  LQR  weights: 

-  (ffu)*''* 

(61) 

(62) 

As  was  shown  earlier,  requiring  the  matrix 


to  be  positive  definite  satisfies  the  necessary  constraints 
for  Ricatti  equation  (56)  to  provide  a  stabilizing  solution. 
In  this  example,  positive  definiteness  of  the  above  matrix 
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provides  the  following  constraints: 


<Jix  >  0 
Ob  >  0 
Qxi  > 

These  constraints  and  equation  (61)  show  that  ki  can  be 
selected  arbitrarily  such  that  ki  >  0.  Also  from  the 
constraints,  Si  can  be  selected  such  that 

-  ^x 

approaches  zero;  therefore,  from  equation  (62),  k,  is  solely 
dependent  on  Om  and  can  also  be  selected  arbitrarily  to 
satisfy  kt  >  0.  As  stated  before,  this  allows  for  assigning 
any  closed  loop  eigenvalues  In  the  left  half  of  the  complex 
plane. 

Recall  that  the  regulator  gains  and  Ricatti  equations 
given  above  are  extensions  of  standard  LQR  theory  and  are 
valid  without  cross  coupling  weights.  To  revert  back  to  the 
standard  LQR  relationships,  Ss[0]  Is  substituted  Into  the 
appropriate  equations.  For  this  example,  forcing  S=[0] 
yields  the  following: 

iCa  -  [fl&a  ♦  (63) 

Notice  from  equation  (61)  that  k,  can  still  be  arbitrarily 
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selected  via  qn,  but  kt  is  now  restricted  by  the  selection 
of  ki.  In  fact,  k,  can  never  be  less  than  (2ki)‘  and  hence 
the  LQR  achievable  region  is  restricted  to  a  closed  loop 
damping  factor  of  greater  than  0.707.  Thus  it  is  easily 
seen  that  the  addition  of  the  cross  coupling  weight,  s,  in 
this  two  state  example,  directly  allows  for  the  arbitrary 
placement  of  the  closed  loop  eigenvalues  within  the  left 
half  complex  plane.  Therefore,  the  restrictions  imposed  by 
the  standard  LQR  solution,  without  cross  coupling,  are 
removed.  Unfortunately,  the  minimum  gain  and  phase  margins 
of  equations  (49)  and  (50)  are  now  no  longer  guaranteed. 
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1}L.  Batulta 


Inputs 

within  this  thesis,  nunerical  data  has  been  rounded  off 
for  ease  of  presentation,  with  a  goal  of  accuracy  to  the 
fourth  significant  digit.  Input  data  must  be  accurate  to 
the  eighth  significant  digit  In  order  to  precisely  duplicate 
the  results  presented  here.  Differences  between  results  In 
this  thesis  with  those  In  Huckabone's  work  [2]  using  the  AH- 
64  model  are  directly  attributable  to  a  difference  in 
accuracy  beyond  the  fourth  significant  digit.  The  AH-64 
model,  represented  In  appendix  A,  Is  displayed  exactly  as 
Input  to  the  algorithm  In  this  thesis.  Inputs  to  the 
algorithm  Include: 


A,  B  -  matrices  from  models 

Ed  -  diagonal  matrix  containing  the  desired 
eigenvalues 

Vd  -  modal  matrix  of  desired  eigenvectors 

(columns  correspond  to  those  of  Ed) 

Fe  -  row  vector  of  weights  corresponding  to 
the  diagonal  elements  of  Ed 

Fv  -  matrix  of  weights  corresponding  to  vd 

tol  -  convergence  tolerance 

r-code  -  1  ^  RM*RM  =  pi 

2  ^  RN*fm  Is  diagonal 

3  ^  RM*RM  Is  general 

s-code  -  0  ^  S  =  [0] 

1  -»  S  Is  filled 


kmax  -  maximum  optimization  Iterations 
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The  desired  closed  loop  elgenstructure  for  both  models 
is  the  same  as  developed  by  Garrard  and  Llebst  [12]  and  is 
presented  in  Table  I.  The  selected  eigenvalues  and 
eigenvectors  were  based  on  work  by  Hoh  [23] »  which  was  a 
precursor  to  ADS-33  [8].  Unity  values  represent  desired 
state-to-state  coupling.  The  non-unity  elements 
corresponding  to  0  and  0  are  the  inverses  of  the  desired 
roll  and  pitch  eigenvalues,  as  these  states  are  the 
integrals  of  the  respective  rates.  An  X  denotes  an  element 
of  the  eigenvector  where  coupling  is  inevitable.  These 
values  are  allowed  to  float  freely  in  the  algorithm  by 
applying  a  zero  weighting  factor  to  the  corresponding 
element  in  the  matrix  Fv. 


Table  I 

Desired  Closed  Loop  Elgenstructure 


State 

Mode  1 

Long 

Vel 

Lat 

Vel 

Heave 

Pitch 

Yaw 

Roll 

X 

-.801  ± 
.3871 

-.802  1 
.3881 

-1.0 

-2.9 

B 

-3.5 

u 

1 

0 

0 

X 

0 

0 

V 

0 

1 

0 

0 

0 

X 

w 

0 

0 

1 

0 

0 

0 

p 

0 

X 

0 

0 

0 

1 

q 

X 

0 

0 

1 

0 

0 

r 

0 

0 

0 

0 

1 

0 

0 

X 

0 

0 

0 

-.2857 

6 

X 

0 

0 

-.345 

0 
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A  small  problem  with  the  algorithm  was  discovered 
during  the  research  to  determine  stability  margins.  When 
desired  eigenvalues  were  Input  as  unstable,  the  algorithm 
did  not  work  properly.  Specifically,  the  algorithm  yielded 
an  unstable  solution.  It  Is  felt  that  this  result  Is  due  to 
a  fault  within  the  DLQGLIB  subroutine  REG  which  allows  for 
an  unstable  solution.  If  an  unstable  solution  Is  provided 
as  a  possible  achieved  solution  to  be  compared  to  an 
unstable  desired  structure,  the  algorithm  may  select  It  as 
the  solution  that  minimizes  3.  In  this  thesis,  unstable 
desired  structures  were  Input  only  In  an  attempt  to 
determine  If  any  achievable  region  could  be  mapped  out  for 
the  LQR  with  cross  coupling  weights.  As  there  are  few 
Instances  where  an  unstable  solution  Is  desired,  this 
problem  does  not  appear  to  provide  any  major  obstacles  In 
the  effective  use  of  the  algorithm. 

As  was  mentioned  before,  the  order  Is  Important  In 
entering  the  desired  elgenstructure.  Figure  2  shows  the 
desired  eigenvalues  and  achieved  eigenvalues  for  the  AH-64 
where  Q  and  R  are  appropriately  dimensioned  Identity 
matrices,  which  Is  the  first  guess  In  the  algorithm.  Some 
Insight  Is  gained  from  this  In  deciding  the  Input  order. 
Specifically,  the  complex  desired  eigenvalues  must  be  paired 
before  the  real  eigenvalues  In  order  to  ensure  that  they  are 
matched  up  with  complex  achieved  eigenvalues.  Also,  the 
desired  eigenvalue  at  (-3.5,  0)  must  be  paired  last  so  that 
It  will  match  up  with  the  achieved  eigenvalue  at  (-18.1,  0) . 


35 


Figure  2 

AH-64  Eigenvalue  Pairing 


Thus  the  desired  helicopter  elgenstructure  for  this  thesis 
was  always  Input  with  the  most  dominant  eigenvalue  first 
proceeding  to  the  least  dominant,  unless  otherwise  stated. 
The  eigenvalue  with  the  positive  Imaginary  component  Is 
always  Input  first. 
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Unity  weighting  of  the  eigenvalues  refers  to  each 
element  of  Fe  being  set  equal  to  one.  Unity  weighting  of 
the  eigenvectors  refers  to  each  element  of  Fv  being  set 
equal  to  one  except  those  elements  corresponding  to  the 
components  of  Vd  that  are  free  to  float  (see  Table  I). 

These  elements  are  always  set  equal  to  zero. 

The  convergence  tolerance  (tol)  Is  used  In  comparing  3 
for  consecutive  Iterations  to  determine  If  a  good  minimum 
has  been  achieved.  It  was  discovered  that  the  program 
occasionally  reaches  plateaus  where  very  small  Improvements 
to  3  are  made  which  can  be  followed  by  much  larger 
Improvements.  This  is  Illustrated  In  Table  II.  Too  high  of 
a  convergence  tolerance  will  prevent  the  program  from 
achieving  these  large  Improvements.  Too  low  of  a  tolerance 
may  not  allow  the  program  a  natural  stopping  point >  thus  a 
maximum  Iteration  value,  kmax,  must  be  set. 

Unless  otherwise  stated,  a  value  of  10**  Is  used  for 
tol  and  150  Is  used  for  kmax.  While  10'*  may  be  beyond  the 
accuracy  desired  In  measuring  the  performance  Index,  there 
were  cases  where  higher  values  for  tol  were  shown  to  mask  a 
27%  reduction  of  3.  And,  while  a  kmax  value  of  150  would 
appear  to  yield  adequate  opportunity  for  the  program  to 
converge,  cases  were  run  to  this  limit  utilizing  over  eight 
hours  of  central  processing  unit  time,  while  continuing  to 
reduce  3  by  more  than  10**.  It  Is  felt,  however,  that  these 
values  provide  a  good  basis  for  comparison  of  results  using 
different  r-codes  and  s-codes. 
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Table  II 

Sample  3  Convergence 


Iteration 

3 

Improvement 

1 

10.961 

2 

9.523 

1.438 

3 

7.703 

1.820 

4 

7.392 

0.311 

5 

7.134 

0.258 

6 

6.943 

0.191 

7 

6.829 

0.114 

8 

6.766 

0.063 

9 

6.533 

0.233 

10 

6.447 

0.086 

11 

6.405 

0.042 

12 

6.272 

0.133 

13 

6.224 

0.048 

14 

6.153 

0.071 

15 

6.039 

0.114 

16 

5.985 

0.054 

17 

5.943 

0.042 

18 

5.577 

0.366 

19 

5.130 

0.447 

20 

4.662 

0.468 

The  r-codes  and  s-codes  are  specified  for  each  example. 
Recall  that  for  the  case  where  S  is  allowed  to  vary, 
s-code  =  1,  the  R  matrix  can  only  be  restricted  to  positive 
definite. 

Analysis  of  RgSMits 

The  algorithm's  eigenstructure  performance  index,  3, 
presented  in  equation  (42)  can  be  broken  down  into  two 
elements  depicting  the  contributions  of  the  eigenvalue  and 
eigenvector  differences.  Unfortunately,  for  weighting 
values  other  than  unity,  a  true  measure  of  the  distance  from 
the  achieved  to  desired  eigenstructure  will  not  be  reflected 
in  the  performance  index  or  its  parts.  In  addition,  the 
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achieved  and  desired  eigenvalues  nay  not  be  properly  paired 
by  the  algorlthe,  as  was  previously  discussed. 

In  order  to  analyze  the  performance  of  the  algorithm 
solutions,  the  parameter  Is  Introduced  as 

•A  -  g 


where  the  parameters  are  paired  by  mode.  This  gives  a  true 
measure  of  the  closeness  of  the  achieved  solution's 
eigenvalues  to  the  desired  eigenvalues.  Note  that  when 
unity  weighting  of  the  eigenvalues  is  used  and  the  returned 
solution  does  properly  match  eigenvalues,  does  reflect 
the  eigenvalue  contribution  to 

Modal  decoupling  is  analyzed  via  the  eigenvectors.  The 
algorithm  normalizes  all  eigenvectors  in  the  program, 
including  the  output.  The  achieved  eigenvectors  presented 
here  have  been  multiplied  by  the  inverse  of  the  eigenvector 
element  corresponding  to  the  primary  desired  response 
element  of  the  desired  eigenstructure.  These  are  the 
elements  valued  at  unity  in  Table  I.  Thus  the  non-unity 
elements  of  each  eigenvector  may  be  viewed  as  a  distance 
away  from  zero  coupling  of  the  respective  element.* 


*The  free-to-f loat  and  non-unity  elements  of  the  desired 
eigenstructure  are  annotated  in  the  data  presentation  tables 
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AH-64  Results 

Achieving  a  Closer  Solution.  Table  III  shows  a  summary 
of  the  results  for  the  AH-64  model  using  unity  weighting. 


Table  III 

Summary  of  AH-64  Results  (Unity  Weighting) 


Run  # 

r-code 

s-code 

3 

a 

1 

1 

0 

4.240 

1 

2 

2 

0 

3.560 

.89371 

3 

3 

0 

3.333 

.47105 

4 

1 

1 

6.219 

.88463 

5 

2 

1 

11.503 

.48946 

6 

3 

1 

2.847 

.80189 

The  first  three  runs  of  this  example  show  expected  results, 
in  that  reduced  restrictions  on  R  allow  for  more  flexibility 
in  reducing  J.  Note  that  the  addition  of  the  cross  coupling 
weights  in  run  6  results  in  a  17%  decrease  of  3.  Runs  4  and 
5  demonstrate  a  problem  with  the  algorithm.  Recall  that 
when  using  a  filled  S  matrix,  s-code  =  1,  the  r-code 
selected  can  not  directly  effect  the  form  of  the  R  matrix  as 
it  does  with  an  S-[0],  s-code  =  0.  Thus  the  r-codes  only 
effect  the  path  taken  by  the  algorithm  in  finding  an  optimal 
solution.  Additionally,  the  algorithm  does  not  provide  for 
the  possibility  of  following  the  same  path  as  with  s-code  - 
0,  which  in  this  case  would  provide  a  better  solution. 

A  surprising  result  from  run  3  is  relatively  poor 
robustness.  The  IGM  and  I PM  associated  with  run  3  are 
[0.680,  1.890]  and  [-27.2®,  27.2®],  respectively.  These 
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Mrglns  Indicate  that  tha  algorlthn  will  not  always  produce 
an  R  Matrix  that  yields  good  robustness  properties,  even 
without  the  use  of  cross  coupling  weights.  This  problen  Is 
addressed  later  In  the  thesis.  Run  6  provides  acceptable 
IGM  and  IPM  of  [0.555,  5.047]  and  [-47. 3",  47. S^], 
respectively.  Minimum  singular  values  for  runs  1,  3  and  6 
are  presented  In  Figure  3  for  comparison. 


frvquwwy  (rod/Mo) 


Figure  3 

AH-64  Minimum  Singular  Values 


The  achieved  elgenstructure  for  runs  3  and  6  are 


presented  In  Table  IV  and  Table  V,  respectively.  The 
appropriate  gain  matrix,  K,  Is  presented  below  each  table 


Table  IV 

AH-64  Achieved  Closed  Loop  Eigenstructure  (Run  3) 


Notes :  ”  denotes  value  free  to  float 

'  denotes  desired  value  of  -0.345 
"  denotes  desired  value  of  -0.286 

K  =  Columns  1  through  4 

4.4662e-01  2.0166e-01  -l.OSlSe-t^OO  -2.6179e-01 

3.1036e-01  -2.2492e-01  4.9597e-01  -8.6887e-02 

3.7330e-02  1.7525e-02  8.8220e-01  9.8005e-02 

-1.6387e-01  -2.5946e-02  6.3103e-01  1.9216e-02 

Columns  5  through  8 

5.6113e-01  1.2189e>00  -8.2418e-02  5.7191e-02 

-7.0436e-01  -5.8664e-02  -1.7295e-01  -8.3843e-01 

7.3459e-02  2.9645e-01  5.0998e-01  -6.8360e-02 

3.9136e-01  -4.9041e-01  -1.4568e-01  -2.2571e-01 
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Table  V 

AH-64  Achieved  Closed  Loop  Eigenstructure  (Run  6) 


Mode  1 

State 

Long 

Vel 

Lat 

Vel 

Heave 

Pitch 

Yaw 

Roll 

X. 

-0.7311 

0.3341 

-0.8231 

0.3021 

-0.281 

-2.941 

-2.859 

-3.636 

u 

1 

0.0651 

0.0471 

-0.046 

-.107" 

-0.012 

0.034 

V 

-0.0941 

0.0751 

1 

-0.037 

-0.018 

0.114 

0.078" 

w 

0.2201 

0.0161 

0.1081 

0.0321 

1 

-0.059 

-0.057 

-0.007 

p 

0.0391 

0.0231 

0.8071 

0.7551" 

-0.009 

-0.144 

1.019 

1 

q 

-0.9501 

0.8971" 

0.0351 

0.0051 

-0.017 

1 

0.393 

-0.224 

r 

-0.0581 

0.2141 

-0.0531 

0.0631 

-0.049 

0.196 

1 

0.597 

-0.0691 

0.0371 

-1.1531 

0.4861" 

0.048 

0.044 

-0.392 

-0.292" 

e 

1.5341 

0.5101" 

-0.0351 

0.0111 

0.069 

-.343" 

-0.155 

0.053 

Notes ;  ”  denotes  value  free  to  float 

'  denotes  desired  value  of  -0.345 
”  denotes  desired  value  of  -0.286 


K  =  Columns  1  through  4 


-2.3441e-01 
7.4176e-01 
1.0283e-01 
-7. 707 40-02 


5.4556e-01 

-2.3506e-01 

1.1148e-02 

3.3675e-01 


-1.1410e-01 

8.1201e-02 

3.5349e-03 

7.8341e-02 


8.6480e-03 

2.0316e-03 

2.4469e-01 

-2.1966e-02 


Columns  5  through  8 


6.1181e-01 
-6.7483e-01 
-7. 66410-02 
2.3170e-01 


6.4257e-01 

-5.1553e-02 

-1.0873e-01 

1.4014e-01 


3.8626e-01 

-1.4942e-01 

5.2847e-01 

3.1558e-02 


5.4422e-01 

-1.0433e-i‘00 

-6.8685e-02 

1.7497e-01 
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While  the  addition  of  the  cross  coupling  weighting  in 
run  6  does  achieve  a  desired  reduction  of  3,  the  eigenvalue 
portion,  increases  from  0.365  to  0.588.  This  increase 
is  due  to  the  placement  of  the  heave  eigenvalue  in  run  6, 
which  accounts  for  88%  of  J^.  In  fact,  all  other 
eigenvalues  are  closer  to  desired  in  run  6  than  in  run  3. 

The  improvements  in  the  overall  3  lie  in  the  eigenvectors, 
particularly  the  heave  associated  eigenvector.  Coupling  is 
tremendous  in  run  3;  in  fact,  vertical  velocity  (w)  is  not 
the  primary  response  element,  as  is  desired. 

As  mentioned  earlier,  a  precise  comparison  between  the 
above  results  and  Huckabone's  [2]  results  for  the  same  model 
was  not  possible;  however,  some  important  differences  were 
discovered.  Of  most  importance  for  this  thesis  is  the  fact 
that  the  addition  of  cross  coupling  weights  did  allow  for  a 
solution  closer  to  the  desired  eigenstructure  than  any 
previous  solutions  without  cross  coupling  weights,  where 
closeness  is  measured  via  3.  Also  notable  was  the  reduced 
stability  robustness  obtained  without  the  addition  of  cross 
coupling  weights  in  run  3. 

Convergence  Tolerance  Effects  on  Stability  Margins. 

The  poor  stability  margins  for  run  3  were  somewhat 
surprising  in  that  they  were  much  worse  than  those 
previously  obtained  in  Huckabone's  work  [2].  While  it  has 
previously  been  shown  that  the  addition  of  cross  coupling 
weights  to  the  algorithm  may  yield  poor  robustness,  run  3 
did  not  Include  this  addition.  Further  research  revealed 
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that  the  convarganca  tolaranca  Input  to  tha  prograa  was 
rasponsibla  for  ralaxad  stability  aarglns. 

Racall  that  tha  LQR  nathod  of  control  systan  dasign  Is 
utilized  primarily  to  take  advantage  of  good  robustness 
properties.  These  properties  were  shown  to  be  dependent  on 
the  type  of  R  matrix  used  within  LQR  theory.  As  stated 
before,  an  R  matrix  close  to  diagonal  could  be  expected  to 
produce  good  stability  margins.  The  algorithm  begins  its 
search  for  an  optimal  solution  by  perturbing  away  from  R=I, 
which  provides  the  guaranteed  stability  margins  presented  in 
equations  (49)  and  (50).  As  R  is  perturbed  away  from  I,  or 
more  correctly  away  from  pi,  R  becomes  far  from  diagonal  and 
the  stability  margins  become  worse. 

In  searching  for  the  lowest  possible  3,  the  convergence 
tolerance  value  was  decreased  significantly  below  the  values 
used  by  Huckabone  [2].  As  mentioned  above,  this  was  done  to 
take  advantage  of  the  apparently  unending  improvements  that 
were  gained  with  the  introduction  of  cross  coupling 
weighting  and  to  provide  for  fair  comparison  of  results. 

What  it  also  did,  was  allow  R  to  perturb  much  farther  away 
from  diagonal  than  it  had  been  previously  allowed,  thus 
revealing  poor  stability  robustness  characteristics. 

In  order  to  demonstrate  this  problem,  run  3  was 
repeated  with  an  increased  tol  value  of  10'*  versus  the 
previous  10'*.  The  solution  yielded  3  »  4.67  and 
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a  =  0.9774.  3  Is  slightly  higher,  but  a  significant 

improvement  is  seen  versus  the  previous  a  =  0.4715.  The  R 
matrix  of  run  3,  where  tol  is  10'*,  is  returned  as 

2.0325  0.4111  0.5768  0.9062 

0.4111  2.6846  -0.0203  -0.2034 

0.5768  -0.0203  3.3425  0.7791 

0.9062-0.2034  0.7791  3.3992. 

The  R  matrix  where  tol  is  10'*  is 

2.0971  0.6235  0.3868  0.4499 

0.6235  2.2980  0.0294  0.0437 
0.3868  0.0294  2.6814  0.6522 
0.4499  0.0437  0.6522  2.1886. 

Notice  that  for  the  lower  tolerance  run,  R  has  a  larger 
spread  between  diagonal  elements  as  well  as  larger  off 
diagonal  elements,  which  is  what  is  meant  by  being  farther 
away  from  diagonal. 

The  loss  of  stability  margins  can  be  even  more 
pronounced  when  cross  coupling  weighting  is  introduced. 
While  the  above  logic  of  perturbing  R  away  from  I  still 
holds,  recall  that  R  can  now  be  perturbed  by  twice  as 
many  variables,  and  may  thus  be  perturbed  away  twice  as 
fast.  This  is  seen  by  reviewing  the  equation 

M  -  tOfMM  *  (65) 
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Nevertheless,  the  stability  Margins  shown  above  for  run 
6  are  not  completely  Inadequate.  In  fact,  stability  Is  not 


a  requirement  for  the  hovering  helicopter  applications 
presented  In  this  thesis  according  to  current  handling 
qualities  requirements  [8:20].  And,  when  run  6  Is  repeated 
with  a  higher  convergence  tolerance  of  10'*,  a  Is  Improved 
from  0.8019  to  0.8612  with  a  resulting  3  -  3.204.  This  3  Is 
still  lower  than  any  previous  runs  for  this  model  without 
cross  coupling  weighting. 

WH-6Q  Rgsgits 

Dimensional  Effects.  The  UH-60  model  was  Input  to  the 
algorithm  with  unity  weighting  on  the  elgenstructure.  A 
summary  of  results  Is  presented  in  Table  VI.  The  algorithm 
was  unable  to  provide  close  matches  to  the  desired 
elgenstructure.  Also,  introduction  of  the  cross  coupling 
weighting  did  not  Improve  the  algorithm's  performance. 

Recall  that,  unlike  the  AH-64  model,  the  UH-60  model  is 
dimensional.  This  prevents  the  algorithm  from  providing 
adequate  results. 


Table  VI 

Summary  of  UH-60  Results  (Unity  Weighting) 


Run  # 

r-code 

s-code 

3 

7 

1 

0 

12.876 

8 

2 

0 

6.693 

9 

3 

0 

7.483 

10 

1 

1 

8.113 

11 

2 

1 

19.707 

12 

3 

1 

9.159 

Because  of  the  absence  of  dimensional  conditioning  of 
the  system,  the  differences  in  units  of  measure  forces  the 
algorithm  to  select  weighting  values  that  make  up  for  the 
differences  in  dimensions.  This  can  be  seen  by  reviewing 
the  resulting  R  matrix  for  the  best  solution,  run  8: 


8.040  000 

0  0 . 826  0  0 

0  0  2.448  0 

0  0  00. 056 


As  would  be  expected  from  the  above  matrix,  robustness  for 
this  solution  is  poor  as  demonstrated  by  a  =  0.523.  The 
requirement  to  vary  the  weighting  matrices  to  account  for 
dimensions  effectively  constrains  the  path  that  the 
algorithm  must  take  in  finding  an  optimal  solution.  Thus 
the  algorithm  reaches  a  dead  end  that  it  would  normally 
avoid  without  the  dimensionally  imposed  constraints.  To 
avoid  this  problem,  the  UH-60  model  was  non-dimensional i zed. 
Changing  scales  of  dimensions  (i.e.  ft/sec  to  m/sec)  could 
be  used  to  provide  the  same  affect. 

The  following  maximum  values  were  used  for  the  given 
states  and  inputs  in  non-dimensional izing  the  UH-60  model: 

60  ft/sec  -  u,  V,  w 

60  ®/sec  -  p,  q,  r 

45  ®  e 

8  inches  -  0,,  0^*  0. 

4  Inches  -  0,. 
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,1,  I.  .MII..U.  I  U.I  .1.1  . . .  .1,1.111.11,1.1111  I.  ...  . . . .  ~ 

For  tho  non>d1iMns1ona1  UH>60,  with  unity  walghting,  tho 
algorlthn  raturnad  a  solution  with  3  »  4.145  and  a  s  0.765 
with  r-coda  »  3  and  s-coda  s  1.  Elganvalua  Hatching  was 
vary  good  with  =  0.185.  Tha  achlavad  alganvactor 
assoc latad  with  tha  haava  moda  Is 


0.2014 
“1.7830 
1.0000 
-3.2322 
“0.4236 
0 . 0643 
3.2208 
0.4188 


This  is  obviously  a  poor  match  accounting  for  ovar  40%  of 
tha  alganvactor  contribution  to  3.  Using  thasa  rasults  as  a 
basal Ina,  savaral  variations  wara  run  to  lllustrata  soma 
dasignar  options. 

Elaenstructura  Input  Ordar  Ef facts.  First,  a  unity 
weighting  casa  for  tha  non-dimanslonal  UH-60  model  was  run 
with  a  different  desired  elgenstructure  order  to  demonstrate 
the  affect  on  the  optimization  path.  The  new  eigenvalue 
ordar  was  arbitrarily  selected  as: 


1. 

“0.801 

- 

0.3871 

2. 

“0.801 

+ 

0.3871 

3. 

“0.802 

- 

0.3881 

4. 

-0.802 

+ 

0.3881 

5. 

“1.0 

6. 

“2.9 

7. 

“3.5 

8. 

“3.0 
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This  run  yielded  the  following  results:  3  =  3.575, 

Jx  -  0.167  and  a  =  0.754.  This  represents  a  14%  decrease  In 
3  with  only  a  slight  reduction  In  stability  margins. 
Obviously,  as  the  program  approaches  a  solution  where  the 
achieved  eigenvalues  are  close  to  the  desired  values,  as  Is 
the  case  for  this  and  the  baseline  runs,  the  eigenvalue 
pairing  will  be  correct.  But,  as  was  shown  before,  the 
Initial  solutions  may  yield  eigenvalues  far  from  desired, 
resulting  In  different  pairings.  This  Initial  pairing 
affects  the  path  that  the  algorithm  follows  In  searching  for 
an  optimal  solution.  For  the  Initial  run  with  the  non- 
dimensional  UH-60  model,  the  solution  was  achieved  after  7 
optimization  Iterations.  The  run  with  the  new  order 
completed  67  Iterations  before  yielding  a  solution.  Note 
that  In  both  cases,  the  standard  convergence  tolerance  of 
10'*  was  used. 

Weighting  Element  Changes.  The  next  change  to  the 
baseline  run  was  a  single  element  weighting  change.  Since, 
the  eigenvector  associated  with  the  heave  mode  was  poorly 
decoupled  In  the  Initial  run,  the  weight  on  the  vertical 
velocity  component  (w)  was  Increased  from  1  to  4.  The 
resulting  solution  yielded 

3  s  2.815,  Jt  s  0.171  and  a  -  0.799.  The  emphasis  on  the 
heave  eigenvector,  through  the  Increased  weight  on  the 
desired  response  element,  allowed  the  algorithm  to  find  a 
closer  solution  with  better  stability  margins.  The  heave 
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•  Igenvactor  was  draiiatically  dacoupled  as  saan  hara 


-0.0353 
-0.0299 
1.0000 
0.0050 
-0.0057 
0.1083 
-0.1228 
0 . 0002 


In  order  to  demonstrate  the  potential  of  the  program, 
several  iterations  of  the  non-dimensional  UH-60  model  were 
run.  The  following  designer  weighting  matrices  were  used: 


wm  •[  2  2  2 


2  3 


2  2  2  ] 


1  1 
1  1 


1.5  1.5 
1  1 
0  0 
1  1 
1  1 
0  0 


1 

1 

3 

0 

1 

2 

0 

1 


1 

1 

3 

0 

1 

2 

0 

1 


1 

1 

6 

1 

1 

1 

1 

1 


0 

1 

2 

3 
1 

4 
1 
1 


1  1 
1  0 
2  1 

4  1 

5  2.5 

6  1 

2  1 
3  1 


The  convergence  tolerance  was  set  at  10‘*  in  an  effort  to 
obtain  good  stability  margins.  The  algorithm  achieved  a 
3  of  2.276.  Note  that  the  algorithm  achieved  a  closer  match 
than  previous  iterations  even  with  non-unity  weights 
included  in  3.  The  minimum  singular  value,  a  =  0.808, 
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yields  an  I PM  of  [47.7®,  -47.7®]  and  an  IGM  of  [0.553, 
5.214].  The  achieved  closed  loop  elgenstructure  Is 
given  In  Table  VII. 


Table  VII 

UH-60  Achieved  Closed  Loop  Elgenstructure 


Mode  1 

State 

Long 

Vel 

Lat 

Vel 

Heave 

Pitch 

Yaw 

Roll 

X 

-0.636± 

0.3051 

-1.0631 

0.4631 

-0.901 

-2.900 

-2.891 

-3.668 

u 

1 

0.096? 

0.0151 

-0.032 

- . 134" 

0.045 

0.020 

V 

0.0631 

0.0391 

1 

-0.029 

0.001 

0.040" 

w 

-0.083T 

0.0291 

0.071? 

0.0891 

1 

0.021 

-0.321 

-0.015 

p 

0.057? 

0.0161 

1.668? 

1.8601" 

0.011 

-0.011 

0.002 

1 

q 

-0.5681 

0.6621" 

-0.1321 

0.1391 

-0.008 

1 

-0.002 

-0.136 

r 

-0.048? 

0.1151 

0.002? 

0.0381 

0.124 

0.009 

1 

0.065 

♦ 

-0.086? 

0.0051 

-1.9621 

0.8981" 

-0.018 

0.004 

-0.017 

-0.274" 

e 

1.134? 

0.4941" 

0.152? 

0.0631 

0.002 

-.345" 

-0.017 

0.036 

Notes ;  ”  denotes  value  free  to  float 

'  denotes  desired  value  of  -0.345 
*  denotes  desired  value  of  -0.286 
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2.  Conclutioni 


The  addition  of  cross  coupling  weighting  to  the 
algorithm  developed  by  Captain  Thomas  Huckabone  [2]  allowed 
greater  flexibility  In  achieving  a  desired  closed  loop 
elgenstructure  for  a  linearized  model  of  a  combat  helicopter 
at  hovering  conditions.  The  designer  Is  able  to  use  the 
algorithm  to  determine  a  full  state  regulator  that  yields 
desired  performance*  decoupling  and  acceptable  robustness. 

Non-dimensional Izat Ion  of  the  Input  system  was  found  to 
be  necessary  In  achieving  good  results  for  the  algorithm. 

The  order  of  the  desired  elgenstructure  was  shown  to  alter 
the  results  of  the  algorithm.  The  convergence  tolerance 
Input  to  the  program  was  shown  to  have  an  Impact  on  the 
stability  margins  of  the  closed  loop  system.  Also*  the 
addition  of  cross  coupling  weights  was  shown  to  remove  the 
guaranteed  stability  margins  of  the  standard  LQR  solution. 
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Recowiiandatlons 


The  algorithm  provides  several  areas  for  further 
research.  These  Include:  program  maintenance,  expanding 
the  application  to  non-hovering  helicopters,  validating  the 
regulator  solution  via  simulation  and  altering  the  algorithm 
to  Improve  Its  robustness  properties.  Some  of  the  following 
proposals  should  be  Incorporated  Into  an  extension  of  this 
and  previous  work  as  a  Master's  thesis  at  AFIT. 

The  source  code  In  appendix  C  has  not  been  properly 
reviewed  with  respect  to  operation  and  can  probably  be  made 
much  more  efficient.  A  good  start  would  be  to  eliminate  the 
use  of  the  DSVR6P  subroutine  from  the  IMSL  package.  In 
addition  to  Improving  efficiency,  this  will  allow  for  the 
program  to  be  used  In  a  stand  alone  capacity  on  any  computer 
system  that  uses  the  FORTRAN  language.  Additionally,  the 
DLQGLIB  package  has  exhibited  some  problems  and  should  be 
altered  or  replaced. 

This  thesis  has  limited  the  application  of  the 
algorithm  to  two  conventional  combat  helicopter  models  In 
hovering  flight.  In  addition  to  looking  at  expanded  flight 
conditions  for  helicopters,  the  algorithm  should  be  applied 
to  other  multi -variable  systems  that  can  benefit  from  the 
application  of  optimal  control  design  techniques. 

Since  full  state  feedback  Is  not  realistic  In  most 
cases,  continued  design  efforts  should  be  carried  out  to 
verify  the  performance  and  stability  capabilities  yielded  by 
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th«  algorlthn  via  tlna  rasponsa  slnulatlons.  Additionally, 
tha  algorithm  could  ba  mada  to  apply  to  obsarvar  dasign  as 
wall  as  regulator  design. 

As  reported,  the  algorithm  currently  does  not  directly 
provide  good  stability  margins  when  applying  cross  coupling 
weights.  The  addition  of  a  stability  parameter  to  the 
algorithm  would  correct  this  problem.  Anderson  and  Moore 
[19: sect.  3.5]  show  that  a  prescribed  degree  of  stability 
can  be  Introduced  to  the  standard  regulator  problem  via  a 
redefined  LQR  performance  Index 

J  "  f  e***  (x'flar  +  a*ko)dt  (66) 

0 

where  o  specifies  the  minimum  degree  of  stability  of  the 
resulting  closed  loop  system.  The  algorithm  could  be 
altered  to  Include  a  as  a  designer  chosen  parameter,  thus 
enforcing  a  robust  algorithm  solution. 
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Appendix  A; _ AH-64  Model 


Matrix  A-: 

Columns  1  through  3 

-2.8600000e>02 

7.7900000e-02 

4.6000000e-03 

5.65837506-01 

3.36637506-01 

2.79337506-01 

0.00000006-t-00 

0.00000006-1-00 

Columns  4  through  6 

3.197207686-03 

-1.157417106-01 

-5.291448576-03 

-2.700000006-^00 

-9.200000006-03 

-1.050000006-1^00 

1.000000006-1-00 

0.000000006-1-00 

Columns  7  through  8 

0.00000006-1-00 

4.46771386-01 

2.23385696-02 

0.00000006-1^00 

0.00000006-1-00 

0.00000006-1-00 

0.00000006-^00 

0.00000006-^00 


-6.37000006-02 

-2.31000006-01 

-2.57000006-02 

-3.58125006-1-00 

8.45175006-01 

-3.50962506-01 

0.000000064^00 

0.00000006-1-00 


1.11274006-01 
-1.43804546-02 
3.14136136-02 
-1.34000006-01 
-7 . 50000006-01 
4.13000006-01 
-5.10000006-03 
9.99000006-01 


2.05000006-02 

5.90000006-03 

-2.61000006-01 

6.80437506-01 

1.43250006-02 

5.73000006-02 

0.00000006-1-00 

0.00000006+00 


-3.58813266-03 

-2.28970336-02 

3.05759166-02 

-6.62000006-01 

2.44000006-02 

-4.00000006-01 

1.03000006-01 

4.99000006-02 


-4.46771386-01 

2.28970336-03 

-4.57940666-02 

0.00000006+00 

0.00000006+00 

0.00000006+00 

0.00000006+00 

0.00000006+00 
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Matrix  B-; 

Columns  1  through  3 


1 . 56600006-01 
-5.69369376-02 
-1.53720006-01 
-1.12940006-t-00 
1.85700006-01 
2.06280006-I-00 
0.00000006-1-00 
0.00000006-»-00 

Column  4 

-6.66666676-04 

2.08700006-01 

8.88888896-04 

4.24020006-1-00 

1.00700006-01 

2.41160006-1-00 

0.00000006-1-00 

0.00000006-1-00 


3.45600006-01 
8.11940306-02 
3 . 45000006-02 
-2.57850006-^00 
-4.34050006-1-00 
4.16900006-01 
0.00000006-(-00 
0 . 00000006-1-00 


-3.99000006-02 

1.71800006-01 

-8.70000006-03 

1.62195006-fOl 

-2.25620006-1-00 

5.01370006^00 

0.00000006-1^00 

0.00000006-1-00 
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PiDDAtldiiL  B; _ UH-6Q  Mathematical  Model 


Matrix  A: 


Columns  1  through  4 


-2.1100e>02 

4.8000e-03 

-5.00000-04 

3.6800e-02 

3.8000e-03 

6.0000e-04 

0 

0 


-1.3200e-02 

-2.0700e-02 

-5.0000e-04 

-3.1200e-02 

6.7000e-03 

4.1000e-03 

0 

0 


1.2300e-02 

0 

-2.3560e-01 

2.3000e-03 

l.lOOOe-03 

0 

0 

0 


-1.2081e-i'00 

-2.5520e-01 

-9.7700e-02 

-5.7728e+00 

1.5320e-01 

-6.3800e-02 

l.OOOOe-i^OO 

0 


Columns  5  through  8 


2.0151e<i'00 

-8.6590e-01 

1.9317e+00 

-1.6579e+00 

-9.0940e-01 

-1.1250e-01 

-2.3000e-03 

9.9870e-01 


Matrix  B; 

5.8590e-01 

1.0980e-01 

-6.4483e-i-00 

-2.1450e-01 

-3.7700e-02 

1.0900e-01 

0 

0 


-2.8160e-01 

-1.3799e+00 

2.2755e+00 

1.4310e-01 

-1.8500e-02 

-2.2380e-01 

4.6200e-02 

5.0500e-02 


-1.4993e-i-00 

-3.7400O-02 

-6.5900e-02 

-9.8000e-03 

3.4830e-01 

1.3200e-02 

0 

0 


-3.6500e-02 
3. 202804^01 
1.6168e-i-00 
4.5500e-02 
1.3000e-02 
l.OOOOe-03 
0 
0 


-7.1600e-02 

4.1700e-01 

-2.2800e-02 

1.3989e-i'00 

2.2800e-02 

2.3300e-02 

0 

0 


-3.2066e-f01 
7.5000e-02 
-1. 482804^00 
0 
0 
0 
0 
0 


9.5370e-01 

-8.7780e-01 

4.3320e-01 

-5.9880e-01 

-9.4000e-03 

3.9980e-01 

0 

0 
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ADD*ndix  ct _ Sgurca  Codt 


PROGRAM  EIGSPACE 
IMPLICIT  C0MPLEX*16  C 
IMPLICIT  REALMS  (A-B»D-H,0-Z) 

COMMON  /INOU/KIN,KOUT 

COMMON  A,B,ed,ea,G,NR,NA,ND,M,N,NN,ACL,P,EV, 

+  WOES, WACH, ca 1 pha, 1 write, nrcode»nscode 
DIMENSION  X(BO) ,A(10,10) ,B(10,10) ,R(10,10) ,Q(10,10) , 
1  RK(10,10) ,G(10,10) ,ACL(10,10) ,P(10,10) ,EV(10) 
DIMENSION  XGUESS ( 80 ) , XS ( 80 ) , GRAD ( 80 ) , X2 ( 80 ) 

DIMENSION  v(80,81) ,Fvec(81) ,vs(80) ,vss(80,81) , 

•f  Fvecl(81) 

DIMENSION  EVS(10),PS(10,10),S(10,10) 

DIMENSION  edr(lO) ,edi (10) ,wdesr(10,10) ,wdes1 (10,10) 
REAL*8  maged(lO) 

C0MPLEX*16  ea ( 10 ) , edg ( 10 ) , WDES (10,10), WACH (10,10) 
C0MPLEX*16  ed(lO) ,W0ESS(10,10) ,ca1pha(10) 

INTEGER  IPERMM(81),IPERMD(10) 

EXTERNAL  PPFUNC,DMACH,DUMCGF 

open (UNIT=10, FI LE=' input.dat' ,STATUS='old' ) 

open ( UNIT=9 , FI LE= ' output . dat  * , STATUS= 'old') 

rewind  10 

rewind  9 

KIN-5 

K0UT=6 


c  read  model  size  and  set  integers  for  array  sizes 


iwrite=0 

readdO,*)  N 

NA=N 

NN=2*N 

NA2sN*N 

N0=NN*(4*N+3) 


c  read  values  for  A  and  8  matrices 


jj=l 
1count-0 
do  10  i-l,NA2 
i  counts  i  count'd 
if (icount.eq.il)  then 
Jj=jj+1 

i count si 
endif 

read (10,*)  A(i count, JJ) 
10  continue 
readdO,*)  M 
NRsM 
NBsN*M 

JJ=1 
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1 countsO 
do  20  i=l,NB 
i counts i count+1 
If (1count.eq.ll)  then 
1 count si 
jj=jj+l 
endif 

readdO,*)  B(icount,jJ) 

20  continue 
c - 

c  read  the  desired  elgenstructure  and  weights 
c - 

do  30  1si,N 

readdO,*)  EV(  1 )  ,edr(  1 )  ,ed1  ( 1 ) 
ed( 1 )sDCMPLX(edr( 1 ) ,ed1 ( 1 ) ) 

30  continue 
jj=l 
1 countsO 
do  35  1s1,NA2 
1  counts  1  count'f  1 
If ( 1count.eq.ll)  then 
1 count si 

JJ=jj+l 

endif 

readdO,*)  P(  Icount,  jj )  ,wdesr(  1count>  jj ) , 

+  wdesi ( Icount) jj ) 

WOES( Icount) jj )sOCNPLX(wdesr( Icount) jj ) ) 

*  wdesi ( Icount) jj ) ) 

35  continue 

call  WNORM(WDES)N) 

c - 

c  read  tolerances  and  codes 

c - 

readdO)*)  tol 
1 eval max s 1000 
read (10)*)  nrcode 
read (10)*)  nscode 
read(lO)*)  kmax 

c - 

c  set  Initial  guess  for  R  and  Q.  Use  the  Identity 
c  matrix  In  both  cases.  Put  the  upper  triangular 
c  portion  of  each  In  XGUESS 

c - 

IxsO 

If (nrcode.eq.l)  then 
XGUESS (1) si. OdO 
Ixsl 
goto  51 
endif 

If (nrcode.eq.2)  then 
do  41  Isi.M 
1xslx+l 

XGUESS(1x)sl.0d0 
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u  u  o 


41  continue 
else 

IcountsO 
do  50  i-l,H 

1  counts  1  count -t-l 
do  40  Jjs1count,M 
ixsix+1 

If ( Icount.eq. jj )then 
XGUESS(1x)si.0d0 
else 

XGUESS(1x)s0.0d0 

endif 

40  continue 

50  continue 
endif 

51  continue 
IcountsO 

do  70  1=1, N 

1  counts  i  count-Hl 
do  60  jj=1 count, N 
1x=1x+l 

If ( icount.eq. jJ )then 
XGUESS(1x)si.0d0 
else 

XGUESS(1x)s0.0d0 

endif 

60  continue 
70  continue 


add  S  matrix  parameters  to  end  of  XGUESS  If  applicable 


If (nscode.ne.O)  then 
do  80  1=1, M 
do  85  jjsl,N 
1x=1x+l 

XGUESS(1x)s0.0d0 
85  continue 
80  continue 
endif 


c  Initialize  X,Fvec,Fvecl,vs  and  vss 


do  91  1=1, lx 
X(1)=0.0d0 
vs( 1 )=0.0d0 

91  continue 
1xpl=1x+l 

do  93  1=l,1xpl 
Fvec(1)=0.0d0 
Fvecl(1)=0.0d0 
do  92  j3=l,1x 
vss(33>1)-O.OdO 

92  continue 
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u  u  u  u  u  u 


93  continue 


c  first  Iteration*  change  XGUESS  to  first  solution  (X) 


CALL  FMINS ( i X  *  XGUESS  *  X , tol , i xpl *  v  *  Fvecl , vs  *  vss  * 
+  IPERMM,ieva1iiax) 
kcountsQ 
do  100  isl,ix 
XGUESS(i)=X(i) 

100  continue 


c  begin  subsequent  iterations 


110  continue 

kcount=kcount'»-l 
do  121  i=l,ix 
X2(i)=0.0d0 
vs( i jsO.OdO 

121  continue 

do  123  1-l,ixpl 
Fvec(i)=0.0d0 
do  122  jj=l,ix 
vss(JJ , i jsO.OdO 

122  continue 

123  continue 

cal  1  FMINS( 1x*XGUESS,X2*to1 * ixpl* v* Fvec*vs, 
■f  vss* IPERHM* ievalmax) 
do  130  i*l*ix 
XGUESS(i)=X2(i) 

130  continue 


check  against  tolerance 


WRITE  kcount*Fvecl(l)*Fvec(l) 

del J=dabs(Fvecl(l)-Fvec(l) ) 

Fvecl (l)=Fvec(l) 

if  (de1J.gt.to1 .and.kcount.1t.50)  goto  110 


final  iteration  using  last  XGUESS 


iwrite=l 

CALL  PPFUNCdx* XGUESS* RJ) 
do  238  i=l*N 

write  (***)  ea(i) 

238  continue 
end 
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SUBROUTINE  FMINS(NX,XGUESS,X,TOL,NXPl,v,Fv«c,vs, 

■f  vss,IPERMM,1«va1iiax) 

IMPLICIT  C0HPLEX*16  C 
IMPLICIT  REAL*8  (A-B,D-H,0>Z) 

DIMENSION  XGUESS(NX) ,X(NX) ,v(NX,NXPl) »Fvac(NXPl) , 

•••  vs(NX) , vss(NX,NXPl)  ,vr(80) » vk(80)  ,ve(80) , 

■i'  vt(80)  ,vc(80)  ,vbar(80) 

INTEGER  IPERMM(NXPl) 

1=100 
icallfsO 
1 count=0 

c - 

c  Build  initial  simplex  near  XGUESS 
c  v( 1 , j )=s1mp1ex  matrix 

c  vs ( 1 ) ^scratch  vector 

c  Fvec( 1 )=funct1on  values  corresponding  to  v(1,j) 

c  col umns 

c - 

xnx=df lotj (NX) 
aa=0 . 5d0 

p=aa* ( dsqrt ( xnx+1 . OdO ) +xnx-l . OdO ) / ( xnx^dsqrt ( 2 . OdO ) ) 
q=aa* ( dsqrt ( xnx+1 . OdO ) -1 . OdO ) / ( xnx^dsqrt ( 2 . OdO ) ) 
do  1010  1=1, NX 
v(1,l)=XGUESS(1) 
vs(1)=v(1,l) 

X(1)=XGUESS(1) 

1010  continue 

leal lf=1cal lf+1 
call  ppfunc(NX,vs,Fv) 

Fvec(l)=Fv 

1=1 

do  1040  Jj=l,NX 
do  1020  kk=l,NX 
vs(kk)=X(kk) 

1020  continue 
1-jJ+l 

do  1030  kk=l,NX 
1f(Jj«eq.kk)  then 
v(kk,1)=vs(kk)+p 
else 

v(kk, 1 )=vs(kk)+q 
endlf 

vs(kk)=v(kk,1) 

1030  continue 

1cal1f=1cal1f+l 
call  PPFUNC(NX,vs,Fv) 

Fvec( 1 )=Fv 
1040  continue 


c -  sort  the  simplex  in  ascending  order 

c  IPERMM(I)  =  vector  of  Index  of  sorted  simplex 

c  sort  Is  In  ascending  order 

c -  vsum(1)  =  summation  of  abs(v(:,1) 


do  1050  1=1,NXP1 
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IPERMH(1)»1 
1050  continue 

ca 1 1  DSVRGP ( NXPl , Fvec , Fvec , I PERMM ) 
do  1070  1  si. NXPl 
do  1060  Jjsl.NX 
vss(jj,1)sv(jj,1) 

1060  continue 
1070  continue 

do  1090  1 si, NXPl 
do  1080  jjsi,NX 

v(  j  j ,  1  )=vss(  j  j ,  IpemimC  1 ) ) 

1080  continue 
1090  continue 
1100  continue 

If ( Icount.gt. levalnax)  goto  1130 
testsO . OdO 
vsumsO.OdO 
do  1120  1s2.NXP1 
do  1110  jjsi,NX 

vsuinsdabs(v(Jj .  1  )-v(Jj>l)  )'fvsuni 
1110  continue 

test sdmax 1 ( test . vsum ) 

1120  continue 

If (test. le.tol )  go  to  1130 


c  Initialize  vr.vk,ve.vt»vs.vss.vc»vbar 


do  1121  Isl.NX 
vr( 1 )s0.0d0 
vkl 1)s0.0d0 
ve ( 1 ) ~0 . OdO 
vt( 1 jsO.OdO 
vs ( 1 ) sQ . OdO 
vc( 1 )s0.0d0 
vbar ( 1 ) sO . OdO 
do  1122  jjsi,NXPl 
vss( 1 .33 )=0.0d0 
1122  continue 
1121  continue 

cal  1  FMINSTEP(v.NX.NXPl.Fvec.vr.vk.ve.vt.vs.vss. 
+  vc.vbar. IPERNM) 

1  counts  i  count'i'l 
goto  1100 
1130  continue 

do  1140  1si,NX 
X(1)=v(1.1) 

1140  continue 
return 
end 
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SUBROUTINE  FNINSTEP( v , NX . NXPl , Fvec , vr, vk , ve , 

4-  vt,v8,vss,vc,vbar, IPERNN) 

IMPLICIT  C0HPLEX*16  C 
IMPLICIT  REAL*8  (A-B,D-H,0-Z) 

DIMENSION  V ( NX » NXPl ) , Fvec ( NXPl ) , vr ( NX ) » vk ( NX ) , 
4^  ve(NX)  »vt(NX)  ,vs(NX)  ,vss(NX,NXPl)  ,vc(NX) , 

4-  vbar(NX) 

INTEGER  IPERMM(NXPl) 
leal  1-0 
a1pha=1.0d0 
beta=0 . 5d0 
ganma=2.0d0 
xnx=df lotj (NX) 
do  2020  1=1, NX 
vb=0.0d0 
do  2010  jj=l,NX 
vb=vb4-v(  1 ,  jj ) 

2010  continue 

vbard  )=vb/xnx 
2020  continue 

do  2030  1=1, NX 

vr(  1  )=vbar(  1 ) 4-alpha* (vbar(  1  )-v(  1  ,NXP1) ) 

2030  continue 

leal  1  =  1  cal  14-1 
call  PPFUNC(NX,vr,fr) 
do  2040  1=1, NX 
vk(1)=vr(1) 

2040  continue 
fk=fr 

If (fr.lt.Fvec(l) )  then 
do  2050  1=1, NX 

ve(  1  )=vbar(  1  )4-gamna*(vr(  1  )-vbar(  1 ) ) 

2050  continue 

1ca11=1ca11  4-  1 
call  PPFUNC(NX,ve,fe) 

If (fe. It. Fvec(l) )  then 
do  2060  1=1, NX 
vk( 1 )=ve(1 ) 

2060  continue 

fk=fe 
end  If 

el  se 

if (fr.ge. Fvec(NXPl) )  then 
do  2070  1=1, NX 
vt(1)=v(1,NXPl) 

2070  continue 

ft=Fvec(NXPl) 

else 

do  2080  1=1, NX 
vt(1)=vr(1) 

2080  continue 

ft=fr 
endlf 
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do  2090  1-1, NX 

vc(  1  )*vbar(  i  )-beta’*'(vbar(  1  )-vt(  i ) ) 
2090  continua 

1  cal  1=1  call -fl 
call  PPFUNC(NX,vc,fc) 
1f(fc.lt.Fvec(NX))  then 
If (fc.ge.fr)  goto  2135 
do  2100  1=1, NX 
vk( 1 )=vc( 1 ) 

2100  continue 
fk=fc 
else 

do  2120  1=2, NX 
do  2110  Jj=l,NX 

v(jj,1)=(v(jj,l)+v(jj,1))/2.0d0 
vs(3j)=v(jj,1 ) 

2110  continue 

leal l=1call+l 

call  PPFUNC(NX,vs,Fv) 

Fvec( 1 )=Fv 
2120  continue 

do  2130  1=1, NX 

vk(1)=(v(1,l)+v(1,NXPl))/2.0d0 
2130  continue 

leal l=1cal 1+1 

call  PPFUNC(NX,vk,fk) 

2135  endlf 
end  If 

do  2140  1=1, NX 
v(1,NXPl)=vk(1) 

2140  continue 

Fvec(NXPl)=fk 
do  2150  111=1,NXP1 
IPERMM(111)=111 
2150  continue 

ca 1 1  DSVRGP ( NXPl , Fvec , Fvec , I PERMM ) 
do  2170  1=1, NXPl 
do  2160  jj=l,NX 
vss(jj,1)=v(jj,1) 

2160  continue 
2170  continue 

do  2190  1=1, NXPl 
do  2180  JJ=1,NX 

v(jj,1)=vss(jJ,IPERMM(1)) 

2180  continue 
2190  continue 
return 
end 
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SUBROUTINE  PPFUNC(NX,X,RJ) 

IMPLICIT  COMPLEX*16  C 
IMPLICIT  REAL*8  (A-B,D-H,0-Z) 

COMMON  /INOU/KIN,KOUT 

COMMON  A»B,ecl,ea,G,NR,NA,ND,M,N,NN,ACL«P,EV, 

WOES ,  WACH ,  cal  pha ,  1  wr  1  te » nrcode ,  nscode 
DIMENSION  X(80),A(10,10),B(10,10),R(10,10),Q(10,10), 

1  RK(10,10),G(10,10)»ACL(10,10),P(10,10),EV(10),evs(10) 
DIMENSION  RCOPY ( 10 , 10 ) , I PVT ( 10 ) , WORK ( 10 ) 

DIMENSION  RIST ( 10 , 10 ) , BRIST (10 , 10 ) , SRIST ( 10 , 10 ) 
DIMENSION  GK(10,10) 

DIMENSION  ANEW( 10, 10) ,QNEW( 10, 10) , BK( 10, 10) 

DIMENSION  0UM(860,1),IDUM(20),WR(10),WI(10),Z(10,10), 

1  IV1(10),FV1(10),ACLS{10,10),S(10,10) 

C0MPLEX*16  ea ( 10 ) , WOES (10,10), WACH (10,10), 

1  WDESS (10,10), WACHS ( 10 , 10 ) , ca 1 pha ( 10 ) 

C0MPLEX*16  ed ( 10 ) , cedi f ( 10 ) , edtmp ( 10 ) , eatmp ( 10 ) 

REAL«8  magea(lO) 


c  this  subroutine  calls  the  cost  function  subroutine, 
c  and  allows  variable  arrays  to  be  set 


RJ=:0.0d0 
do  176  1=1,10 

ea(1)=DCMPLX(O.OdO,O.OdO) 
do  177  jj=l,10 
ACL(1,jj)=0.0d0 
WACH  ( 1 ,  J  j )  =dciiip  1 X  ( 0 .  OdO ,  0 .  OdO ) 

G(i,jj)=0.0d0 
177  continue 
176  continue 

call  PP(NX,X,RJ,NR,NA,ND,N,M,NN,A,B,R,Q,S,RK,G,ACL,P, 

1  EV,evs,OUM,IDUM,WR,WI,Z,IVl,FVl,ea,ed,WDES,WACH, 

2  cedi f , edtmp , eatmp , magea , ACLS , WDESS , WACHS , cal pha , 

3  1 wr 1 te , nrcode , nscode , RCOPY, I PVT , WORK, RIST, BRIST , 

4  SRIST,GK,ANEW,QNEW,BK) 

RETURN 

END 
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SUBROUTINE  PP(NX,X,RJ,NR,NA,NO,N,M,NN, A,B,R,Q,S,RK,G, 

1  ACL,P,EV,evs»OUM,IDUM,WR,WI,Z,IVl,FVl,ea,ed>WDES, 

2  WACH , cedi f , edtmp , eatmp , magea , ACLS , WDESS , WACHS , cal pha , 

3  1 wri te , nrcode » nscode , RCOPY , I PVT , WORK , RIST, BRIST , 

4  SRIST,GK,ANEW,QNEW«BK) 

IMPLICIT  C0MPLEX*16  C 
IMPLICIT  REALMS  ( A-B,D-H,0-Z ) 

COMMON  /INOU/KIN,KOUT 

DIMENSION  X(NX),A(N,N),B(N,M)»R(M,M),Q(N,N), 

1  RK(N,N).G(M,N),ACL(N,N)»P{N,N),EV(N),WN0RMA(10), 

2  PS(10,10),EVS(N),S(N,M) 

DIMENSION  DUM(ND,1) , IDUM(NN) ,WR(N) ,WI (N) ,Z(N,N) > 

1  IV1(N),FV1(N) 

DIMENSION  RM(10,10),QH(10»10)»SN(10,10),SNT(10»10) 
DIMENSION  S1(10,10)«S2(10,10),QS{10,10),SR(10,10) 
DIMENSION  er(10),e1(10) 

DIMENSION  ediffflag(10),ACLS(N»N) 

DIMENSION  PDUM( 10,10), EV0UM( 10,10) 

DIMENSION  RC0PY(M,M) , IPVT(M) ,W0RK(M) 

DIMENSION  RIST(M,N) ,BRIST(N,N) ,SRIST(N,N) 

DIMENSION  GK(M,N) 

DIMENSION  ANEW(N,N) ,QNEW(N,N) ,BK(M,N) 

D I MENS I ON  QSAVE (10,10), RSAVE (10,10) 

C0MPLEX*16  ea(N) ,ed(N) ,WDES(N,N) ,WACH(N,N) ,WDESS(N,N) , 
1  WACHS(N,N),ca1pha(N) 

C0MPLEX*16  cedlf (N) ,edtmp(N) ,eatnp(N) 

INTEGER  IPERMA(10),im1n(10) 

REALMS  magea (10) 

LOGICAL  ELIM 
IPRT=0 


c  if  last  Iteration,  write  to  output 


if ( i write. ne.O)  then 
write  (9,*)  N 
write  (9,*)  M 
endif 


set  the  Q,  R  and  S  matrices 


CALL  MAKEQRS(N,M,NX,X,Q,R,S,RM,QH,SN,SNT,S1,S2,QS,SR, 
1  nrcode, nscode, QSAVE, RSAVE) 
if ( iwrite.ne.O)  then 
do  557  i=l,M 
do  556  jj=l,M 

write  (9,*)  R(jj,i) 

556  continue 

557  continue 

do  555  i-l,N 
do  554  Jj=l,N 

write  (9,*)  Q(jj,i) 

554  continue 

555  continue 
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QUO  U  U  U 


do  553  1sl»M 
do  552 

write  (9,*)  S(jj,1) 

552  continue 

553  continue 
endlf 


c  calculate  the  Iqr  gain  matrix,  GK 


CALL  TRNATB(N,M,N,M,S,RIST) 

CALL  SAVE(M,M,M,M,R,RCOPY) 

CALL  MLI NEQ ( M , M , N , RCOPY , R 1ST , NCOND , I PVT , WORK , 1 ) 

CALL  MMUL(N,M,N,N,M,N,B,RIST,BRIST) 

CALL  MMUL(N,M,N,N,M,N,S,RIST.SRIST) 

CALL  MSUB(N,N,N,N,N,A,BRIST,ANEW) 

CALL  MSUB(N,N,N,N,N,Q,SRIST,QNEW) 

CALL  REG(NA,NR,N,M,NN,ANEW,B,QNEW,R,RK,G,ACL,DUM,IDUM, 
1  IPRT) 

CALL  MAOD(M,M,M,M,N,G,RIST,GK) 

If ( 1 write. ne.O)  then 
do  530  1=1, N 
do  5555  jj=l,M 

write  (9,*)  GK(jj,1) 

5555  continue 
530  continue 
endlf 


calculate  the  new  closed  loop  eigenvalues 


CALL  MMUL(N,M,N,N,M,N,B,GK,BK) 

CALL  MSUB(N,N,N,N,N,A,BK,ACLS) 
ipc=l 

CALL  EIGVV(NA,N,ACLS,WR,WI,Z,IV1,FV1,IPC,IERR) 


count  complex  pairs  and  configure  eigenvectors 


1comp1ex=0 
do  539  1=1, N 

1f(WI(1 ) .ne.O.OdO)  1 comp! ex=i comp! ex+1 
539  continue 

NCMP=N-1 complex/2 
11=0 

do  540  1=1,NCHP 
11=11+1 

If (dabs(WI ( 1 1 ) ) .gt.0.0)  then 
do  535  jj=l,N 

WACH(jj,11)=DCMPLX(Z(jj,11),Z(jj,i1+l)) 
WACH( j j , 1 1+1) =0C0NJG(WACH( j j ,11)) 

535  continue 

11=11+1 
else 

do  536  jj=l,N 

WACH ( j  j , 1 1 ) =DCMPLX ( Z ( j  j , i 1 ) , 0 . OdO ) 
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536  continue 

end  if 

540  continue 

call  WN0RM(WACH,N) 
do  30  1:^1, N 

ea ( 1 ) =dcmp1 x ( WR ( 1 ) » WI ( 1 ) ) 

30  continue 
c - 

c  write  to  output  if  final  Iteration 
c - 

40  continue 

If ( IwrIte.ne.O)  then 
do  43  1=1, N 

earea1=drea1 (ea(  1 ) ) 
ea1mag=d1mag(ea(  1 ) ) 
write  (9,*)  eareal 
write  (9,*)  ealmag 
43  continue 
endlf 

do  41  1=1,  N 
eatmp( 1 )=ea( 1 ) 
edtmp ( 1 ) =ed ( 1 ) 

EVS(i)=EV(1) 
do  42  jj=l,N 

WACHS(jj,i)=WACH(jj,i) 

WDESS(jj,i)=WDES(jj,1) 

PS(jj,i)=P(jj,1) 

If ( IwrIte.ne.O)  then 
wrea1=drea1 (wach( jj , 1 ) ) 
w1mag=d1mag(wach( j j , 1 ) ) 
write  (9,*)  wreal 
write  (9,*)  wimag 
endlf 

42  continue 

41  continue 
jjj=0 

do  501  jj=l,N 
do  502  kk=l,N 

If (P( jj ikk) .ne.O.OdO)  jjj=jjj+l 
502  continue 
501  continue 
NL=N 
RJ=0.0d0 
do  50  1=1, N 

c - 

c  calculate  difference  between  desired  eigenvalues  and 
c  achievable  eigenvalues,  algorithm  matches  the  closest 

c  achieved  eigenvalue  to  the  desired  eigenvalues.  In  the 

c  order  they  were  Input, 
c - 

do  51  jj=l,NL 

cedlf ( jj )=edtmp(l)-eatmp( jj ) 

51  continue 
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do  56  jJ=l,NL 

ed1fraag(3j )=dsqrt(drea1 (cedif ( jj ) )**2* 

+  diraag ( cedif (jj ) )**2) 

56  continue 

call  iroins(NL,edifmag,  1niin(  1 ) ) 
ii1=IMIN(i) 

RJVEC=O.OdO 
if(jjj.eq.O)  then 
if ( iwrite.ne.O)  then 
goto  503 
end  if 
endif 

c - 

c  calculate  the  calpha  for  the  eigenvectors 
c - 

suml=0.0d0 

sum2=0.0d0 

if (dimag(eatmp( i i i ) ) .ne.O.OdO)  then 
do  584  jj=l,N 

suml  =  suml-t-dimag(WDESS(  jj ,  1) ) 

+  *dreal (WACHS(jj,iii)) 

+  -dreal (WDESS( j j , 1) )*dimag(WACHS( j j , i i i ) ) 

suin2  =  sum2-t-drea1  (WDESS(  jj  ,1) ) 

+  *dreal (WACHS( jj , i i i ) ) 

+  +dimag(WDESS( j j , 1) )*dimag(WACHS( j j , i i i ) ) 

584  continue 
phl3:datan2  ( suml ,  sum2 ) 
ca1phalsdcmp1x(dcos(phl) ,dsin(phl) ) 
cal pha2s-l*ca1 phal 

else 

ca1 phal- ( 1 . OdO , 0 . OdO ) 
ca1pha2=-l*calphal 
endif 

c - 

c  determine  which  calpha  produces  minimum  cost  and 
c  calculate  eigenvector  contribution  to  performance 

c  index 

c - 

DELWIl=0.0d0 
DELWI2=0.0d0 
do  585  jj=l,N 

DELWIl=PS(jj,l)*((DREAL(WDESS(jj,l)- 
+  ca1phal*WACHS(jj,iii)))**2 

+  +(DIMAG(WDESS(jj,l)-calphal 

+  *WACHS(jj,iii)))*’*'2) 

+  +DELWI1 

DE  LW 1 2= PS ( j  j , 1 ) * ( ( DREAL ( WOESS ( j  j , 1 ) - 
+  calpha2*WACHS(jj,iii)))**2 

+  •l•(DIMAG(WDESS(jj,l)-ca1pha2 

+  *WACHS(jj,iii)))*»2) 

+  +DELWI2 

585  continue 
if(DELWIl.lt.DcLWI2)  then 
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DELWI-DELWIl 
cal phal seal phal 
else 

DELWI=DELWI2 
cal phai seal pha2 
endif 

RJVEC=DELWI 
If  ( iwrIte.ne.O)  then 
do  587  jj=l,N 

if (eatmp( Hi) .eq.ea( jj ) )  ca1pha( j j )sca1phai 
587  continue 

if ( 1 .eq.N)then 
do  586  jJsl.N 

alrealsdreal (ca1pha( jj ) ) 
a1 imag=dimag(ca1pha( j j ) ) 

586  continue 

endif 
endif 

503  continue 


c  add  eigenvalue/eigenvector  contribution  to  performance 

c  index 


RJ=RJ+EVS ( 1 ) *edi f mag ( i i i ) **2+RJVEC 


c  reset  the  eigenvalue/vector  arrays  to  eliminate 
c  those  poles  already  matched. 


ksO 

NL=NL-1 
do  52  jj=l,NL 
k=k+l 

ELIMsjj.eq.iii 
IF(ELIM)  K=k+1 
eatmp( jj )=eatmp(k) 
edtmp( j j ) =edtmp( j j+1) 
EVS{jj)=EVS(jj+l) 
if(jjj.eq.O)  goto  591 
do  590  kksl,N 

WACHS ( kk  J  j ) sWACHS ( kk , k ) 
WOESS  ( kk  >  J  J  )  sWOESS  ( kk ,  j  j  -t- 1 ) 
PS(kk,jj)=PS(kk,jj+l) 

590  continue 

591  continue 
52  continue 
50  continue 

RETURN 

END 
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SUBROUTINE  IMINS(NL,EOIFNAG, I ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  EDIFMAG(NL) 

do  5000  jj=l,nl 

If  (dabs(ed1froag(  jj  ) ) .  lt.dabs(edifniag(  i ) ) )  i=jj 
5000  continue 
return 
end 
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SUBROUTINE  WNORH(WVEC,N) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  WNORMV(IO) 

COMPLEX*16  WVEC(N,N) 
do  5  1^1, N 
WNORMV(1)=O.OdO 
5  continue 
do  10  1-1, N 
do  20  jj=l,N 

WNORMV  ( 1 )  ^WNORMV  ( 1 )  -t-  ( drea  1  ( WVEC  ( j  j  ,  1 ) ) )  **2+ 
+  (d1mag(WVEC{jj,1)))**2 

20  continue 
10  continue 
do  30  1-1, N 
do  40  jj=l,N 

WVEC(  jj ,  1  )=WVEC(  jj  ,  1  )/dsqrt(wnorniv(  1 ) ) 

40  continue 
30  continue 
return 
end 
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SUBROUTINE  HAKEQRS(N,M,NX,X,Q,R,S,RM,QH»SN,SNT,S1,S2, 
1  QS « SR , nrcode , nscode , QSAVE , RSAVE ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  X(NX) )Q(N,N) »R(M,M) ,S(N»M) 

DIMENSION  RM(M,M) ,QH(N,N) ,SN(N,M) ,SNT(M,N) 

DIMENSION  S1(N,N) ,S2(M,M) ,QS(N,M) ,SR(N,M) »QSAVE(N,N) , 
DIMENSION  RSAVE(M,M) 

c - 

c  this  subroutine  reads  the  upper  triangular  portion 
c  of  matrices  RM,  QH  and  all  of  SN,  If  applicable,  from 

c  X  and  returns  Q,  R  and  S  matrices  that  meet  the 

c  constraints  necessary  for  an  LQR  solution 
c - 

1x=0 

1 f ( nrcode . eq . 1 . or . nrcode . eq . 2 )  then 
1x=l 

do  116  1=1, M 
do  117  jj=l,M 
If(l.eq.jj)  then 
RM(1,jj)=X(1x) 
else 

RM(1,jj)=0.0d0 
end  If 

117  continue 

If (nrcode. eq.2)  1x=1x+l 
116  continue 
else 

1count=0 
do  101  1=1, M 
1  counts 1 count+1 
do  102  jj=1 count, M 
1x=1x+l 

RM(1,jj)=X(1x) 

RM(jj,1)=X(1x) 

102  continue 
101  continue 
end  If 
1count=0 
do  111  1=1, N 
1  counts  Icount-M 
do  112  jj=1 count, N 
1x=1x+l 

0H(1,jj)=X(1x) 

QH(jj,1 )=X(1x) 

112  continue 
111  continue 

If (nscode.ne.O)  then 
do  131  1=1, M 
do  132  jj=l,N 
1x=1x+l 

SN(jj,1)=X(1x) 

SNT(1,jj)=X(1x) 

132  continue 
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131  continue 
endif 

call  MMUL(N«N,N,N,N,N,QH,QH,Q) 
call  MMUL(M,M,M,M,M,M,RM,RM,R) 

If (nscode.eq.O)  then 
call  ZPART(N,N,M,S) 
else 

call  MHUL(N,M,N,N,M,N,SN,SNT,S1) 
call  MMUL(M,N,M,M,N,M,SNT,SN,S2) 
call  MMUL(N,N,N,N,N,M,QH,SN,QS) 
call  MMUL(N,M,M,N,M,H,SN,RH,SR) 
call  SAVE(N,N,N,N,Q,QSAVE) 
call  MADD(N,N,N,N,N,QSAVE,S1,Q) 
call  SAVE(M,M,M,M,R,RSAVE} 
call  MA0D(M,H,M,M,M,S2,RSAVE,R) 
call  MADD(N,N,N,N,M,QS,SR,S) 
endif 
return 
end 
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ABPgndlX  D; _ OPTatlna  Instructions 


General 

The  following  Is  provided  as  a  guide  for  using  the 
included  algorithm  at  AFIT.  While  the  source  code  is 
written  in  a  standard  FORTRAN  language  [24],  the  program 
requires  some  special  handling.  Particularly,  the  specific 
subroutines  called  upon  in  the  program  must  be  properly 
accessed.  The  procedures  presented  here  ensure  proper  use 
of  the  program.  The  program  is  run  on  the  Virtual  Memory 
System  (VMS)  cluster  at  AFIT.  The  cluster  includes  three 
host  computers:  Hercules,  Lancer  and  Sabre.  Reference  [25] 
contains  more  details  about  the  AFIT  computer  services. 
CgffiPllinfl 

The  EIGSPACE  source  code  is  compiled  using  the  FORTRAN 
compiler  on  VMS.  The  object  file  is  linked  to  the  0LQGLI6 
and  IMSL  subroutine  packages.  The  DLQGLIB  source  code  is 
available  at  AFIT.  For  this  project,  the  DLQGLIB  object 
file  was  placed  in  the  same  directory  as  the  EIGSPACE 
program.  The  IMSL  subroutines  must  also  be  linked  to  the 
algorithm.  These  routines  are  not  available  in  FORTRAN  at 
AFIT.  Following  are  the  commands  used  to  compile  the 
algorithm: 

FORTRAN  eigspace 

LINK  eigspace,dlqg1 ib, imsl ib_share/opt 
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MATLAB""  Interface 


The  EIGSPACE  program  requires  a  single  input  file  and  a 
single  output  file.  Since  most  of  the  input  and  output  are 
expressed  as  matrices »  MATLAB"*  provides  a  convenient 
interface  for  the  program.  MATLAB”*  can  be  used  as  a  stand 
alone  interface,  via  the  LQREA.M  file,  which  requires 
interactive  use  of  a  terminal.  The  program  can  be  more 
efficiently  run  when  submitted  via  the  following  command 
file  EIGSPACE.COM: 


$run  EIGSPACE 
$exit 

The  file,  EIGSPACE.LOG  is  created  which  contains  all  screen 
output  data  from  the  algorithm  as  well  as  other  information 
about  the  program  execution.  The  LQRSAVE.M  file  is  used  to 
create  the  input  file  for  EIGSPACE,  while  the  LQROUT.N  file 
is  used  to  reformat  the  output  file.  The  input  and  output 
files  should  be  in  the  same  directory  as  the  EIGSPACE  file. 
The  m>fi1es  presented  below  ensure  that  the  program  will  run 
properly  when  MATLAB^"  is  used  in  the  same  directory. 
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LQREA  irfllg 


f unction [Q, R, S, ea, va, K] sLQREA(a, b, ed, Fe, vd, Fv, tol ,rcode,scode, 

kmax) 

%LQREA  Elgenstructure  assignment  using  the  Linear  Quadratic 
%Regu1ator. 

% 

%Form  Is  [Q,R,S,ea,va»K]= 

%  LQREA(a,b,ed«  Fe>vd,Fv«to1 , rcode , scode , kmax ) 

% 

% Input  parameters, 

%  a,b=state  space  matrices  of  a  linear  system  with  n  states 
%  and  m  controls 

%  ed=nxn  diagonal  matrix  of  desired  eigenvalues 
%  Fe=nxl  matrix  weighting  each  eigenvalue 
%  vd=nxn  matrix  whose  columns  are  the  desired  eigenvectors 
%  (must  be  In  same  order  as  associated  eigenvalues) 

%  Fv=nxn  matrix  whose  elements  weight  corresponding  elements 
%  of  vd 

%  to1 =convergence  tolerance  for  performance  Index 
%  rcode=(l)  R=ro’»fI,  (2)  Rs[d1ag],  (other)  R=pos1t1ve  definite 
%  scode=:(0)  S=nxm  zero  matrix,  (other)  S=contro1 -state 
%  weights 

%  kmax=max1mum  optimization  Iterations 
iOutput  parameters , 

%  ea=nxn  diagonal  matrix  of  achievable  eigenvalues 
%  va-nxn  matrix  of  achievable  eigenvectors 
%  Q,R,Ssf1na1  state  and  control  weighting  matrices 
%  K-gaIn  matrix 
% 

[n,nc]-s1ze(a) ; 

[mr,m]=s1ze(b) ; 
at=a( : ) ; 
bt=b( : ) ; 
for  1=l;n, 

edd(1,l)=Fe(1); 
edd( 1 ,2)=rea1 (ed( 1,1)); 
edd( 1 ,3)=1mag(ed( 1,1)); 
end 

vdd( ; ,l)=Fv( : ) ; 
vdd( ; ,2)=real (vd( ; ) ) ; 
vdd ( : , 3 ) = 1 mag ( vd ( : ) ) ; 

save  1nput.dat  n  at  m  bt  edd  vdd  tol  rcode  scode  kmax  /ascii 
!run  eigspace 
load  output.dat 
count-1; 

nsoutput( count) ; 

count=count+l; 

m=output( count) ; 

count^count-t-l; 

n2=n*n; 

m2=m*m; 

nm=n*m; 
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Rszeros(m,ii) ; 

R  ( : )  =output  ( count :  count-t'm2-l ) ; 
count=count-*‘ffl2 ; 

Qszeros(n»n) ; 

Q( : )  =output (count: count-»-n2-l) ; 
count-count'i'n2 ; 

S=ones ( n , m ) ; 

S  ( : )  ^output  ( count :  count-t-nm-l ) ; 
countscount-fnm ; 

Kszeros(m,n) 

K ( : )  ^output  ( count :  count-t-nm-l ) 
countscount-fnm 
ea=zeros(n,n) ; 
for  i=l:n 

ea  ( i ,  1 )  ^output  ( count )  j  ^output  ( count-i-l ) ; 
count=count'f2: 
end 

va=zeros(ntn) ; 
for  jj=l:n 
for  i=l:n 

va(  i ,  jj  )=output (count )*t*j*output(count'»-l) 
count=count+2; 
end 
end 
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f unct 1 on  LQRSAVE ( a , b , ed , Fe , vd « Fv , tol , rcode » scode , kmax ) 

%LQRSAVE  saves  all  required  data  for  the  Elgenstructure 
iassignment  algorithm  using  the  Linear  Quadratic  Regulator. 

% 

%  Form  Is  LQRSAVE ( a, b, ed, Fe, vd, Fv,to1 , rcode, scode) 

%  Input  parameters, 

%  a,b  =  state  space  matrices  of  a  linear  system  with 

%  n  states  and  m  controls 

%  ed  =  nxn  diagonal  matrix  of  desired  eigenvalues 

%  Fe  =  nxl  matrix  weighting  each  eigenvalue 

%  vd  =  nxn  matrix  whose  columns  are  the  desired 

%  eigenvectors  (must  be  In  same  order  as 

%  associated  eigenvalues) 

%  Fv  =  nxn  matrix  whose  elements  weight 

%  corresponding  elements  of  vd 

%  tol  s  convergence  tolerance  for  performance  Index 

%  rcode  =  (1)  R  =  pi 

%  (2)  R  =  Idlag] 

%  (3)  R  =  positive  definite 

%  scode  =  (0)  S=[0] 

%  (1)  S=contro1 -state  weights 

%  kmax  =  maximum  optimization  Iterations 

[n,nc]ss1ze(a) ; 

[mr,m]=s1ze(b) ; 
at=a( : ) ; 
btsb( : ) ; 
for  i=l:n, 

edd(1,l)sFe(1); 
edd( 1 ,2)srea1 (ed( 1,1)); 
edd(  1 ,3)siiiiag(ed(  1,1)); 
end 


vdd( : ,l)=Fv( : ) ; 
vdd( : ,2)=rea1 (vd( : ) ) ; 
vdd( : ,3)=1mag(vd( : ) ) ; 

save  1nput.dat  n  at  m  bt  edd  vdd  tol  rcode  scode  kmax  /ascii 


LOROUT  m-file 


function  [Q,R,S,ea,va,K]=LQROUT 

%LQROUT  loads  output  from  EIGNEW  elgenstructure  assignment 
%  algorithm  using  the  Linear  Quadratic  Regulator. 

% 

%  Form  Is  [Q,R,S,ea» va,K]=LQROUT 
% 

%  Output  parameters, 

%  Q,R,S=f1na1  state  and  control  weighting  matrices 

%  easnxn  diagonal  matrix  of  achievable  eigenvalues 

%  va=nxn  matrix  of  achievable  eigenvectors 

%  K=ga1n  matrix 

% 

load  output.dat 
count=l; 

n=output( count) ; 
count=count+l; 
m=output( count) ; 

count^count-fi; 

n2=n*n; 

m2=m*m; 

nm=n*m; 

R=zeros(m,m) ; 

R( :  )=output(  count  :count'fm2-l) ; 
count=count-»-m2 ; 

Qszeros(n,n) ; 

Q( ; )soutput( count :count+n2-l) ; 
countscount*»>n2 ; 

Stones (n,ro) ; 

S( :  )soutput(  count  :count-fnm>l) ; 
count=count+nra ; 

K=zeros(m,n) 

K  ( : )  =output  ( count :  count -t-nm- 1 ) 
count^count+nm 
ea=zeros(n,n) ; 
for  i=l:n 

ea ( 1 , 1 ) =output ( count ) ^output ( count+l ) ; 
countscount't-2; 
end 

va=zeros(n,n) ; 
for  jj=l:n 
for  1=l:n 

va  ( 1 ,  J  j )  ^output  ( count )  ■•■J  ^output  ( count+1 ) ; 
count~count-t'2 ; 
end 
end 
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